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Final  Report 


MODELING  AND  COMPUTATIONAL  ANALYSIS  OF  MULTISCALE 
PHENOMENA  IN  FLUID-STRUCTURE  INTERACTION  PROBLEMS 

1.  Introductory  Comments 

document  summarizes  research  done  during  the  period  June  14,  1989  through 
December  14,  1991  on  a  project  aimed  at  the  development  of  advanced  computational 
methods  for  modeling  multiscale  fluid-structure  interaction  phenomena. 

The  major  thrust  of  the  work  reported  is  the  development  of  the  component  parts  of 
a  new  family  of  computational  schemes  that  could  be  used  to  study  large-scale  problems  in 
fluid-structure  interaction.  The  overall  theme  of  the  project  was  the  optimization  of  the 
computational  process  through  the  development  of  new  adaptive  techniques,  data 
structures,  a  posteriori  error  estimates,  and  high-order  solvers.  The  premise  behind 
studying  these  types  of  high-order  methods  is  that,  with  proper  data  structures,  they  could 
lead  to  exponential  rates  of  convergence  and  thereby  allow  one  to  numerically  solve  very 
large  problems,  with  features  covering  many  different  physical  scales,  using  many  fewer 
degreees  of  freedom  than  that  required  by  conventional  methods.  In  this  spirit,  the  project 
was  an  ambitious  one  for  it  involved  developing  not  only  new  adaptive  algorithms  for 
simulating  the  motion  of  elastic  structures,  but  also  techniques  for  fluid  structure  interaction 
and  for  modeling  various  regimes  of  incompressible  viscous  flows.  Presumably,  these 
methods  could  also  be  used  to  model  the  transition  to  turbulence  and  ultimately  turbulent 

boundary  layers  on  compliant  surfaces.  _ _ _ _ _ 

Particular  approaches  explored  during  this  project  deviated  significantly  from  many 
of  the  linearized  or  quasi-linearized  approaches  prevalent  in  the  signal-processing  literature. 
The  theme  here  was  to  explore  methods  which  attempted  a  direct  assault  on  the  governing 
equations.  While  this  class  of  problems  is  exceedingly  difficult,  the  payoff  could  be 
substantial  in  that  not  only  would  very  general  modeling  capabilities  be  developed,  but  also 
very  large  simulations  could  conceivably  be  handled  by  new  methodologies  which  were 
designed  to  optimize  the  computational  process. 

1.1  State  of  the  Art  at  Start.  The  current  state-of-the-art  in  computational  modeling 
of  large-scale  subsonic  flow-structure  interaction  problems  is  fairly  well-known.  Today  a 
number  of  large-scale  computer  codes  are  available  that  are  typically  based  on  finite 
difference  approaches.  INS  3D,  for  example,  a  NASA— supported  code  developed  for 
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incompressible  flow  simulation  uses  finite  difference  algorithms  and  splitting  techniques  to 
handle  three-dimensional  incompressible  flow  problems  for  Reynolds  numbers  of  around 
1,000.  With  the  addition  of  various  turbulence  models,  Reynolds  numbers  up  to  10^  or 
10^  can  be  handled,  but  the  philosophy  in  resolving  small-scale  flow  features  is  traditional: 
detailed  features  of  the  flow  are  modeled  by  incorporation  of  low  order  schemes  involving 
many  grid  points.  In  complex  flow  calculations,  the  quality  of  a  simulation  is  frequently 
measured  in  terms  of  the  number  of  grid  points  in  the  mesh  and  efficiency  is  guaged  by  the 
number  of  iterations  per  grid  point  per  second.  Few  of  these  calculations  purport  to  model 
detailed  features  of  boundary  layers  or  transition  to  turbulence.  Alternatively,  spectral 
methods  have  been  used  in  recent  years  with  some  success  in  modeling  various  turbulence 
phenomena.  Unfortunately,  these  methods  are  primarily  restricted  to  simple  geometries 
and  periodic  boundary  conditions  and  have  not  been  used  for  significant  interaction 
calculations. 

The  modeling  of  the  structure  itself  and  the  fluid-structure  interface  also  represents  a 
class  of  problems  that  are  at  the  cutting  edge  of  computational  mechanics.  Contemporary 
methods  are  incapable  of  capturing  sufficient  detail  to  adequately  model  features  of 
structure  response  that  influence  such  phenomena  as  acoustics,  or  design  of  sensor 
locations  on  moving  submerged  vehicles,  such  as  submarines,  within  a  turbulent  boundary 
layer.  These  types  of  modeling  probems  are  still  outside  the  reach  of  most  conventional 
methods.  • 

It  has  been  known  for  a  number  of  years  that  in  certain  simple  cases,  such  as  one¬ 
dimensional  elliptic  problems,  a  proper  distribution  of  mesh  parameters  can  produce 
exponential  convergence  of  numerical  simulations;  that  is  to  say,  with  a  decrease  of  the 
mesh  size  or  increase  of  the  spectral  order  or  proper  location  of  grid  points,  the  total 
approximation  error  may  decay  exponentially  as  the  number  of  degrees  of  freedom  are 
increased.  This  fact  led  to  the  question  of  whether  or  not  a  new  family  of  algorithms  could 
be  developed  which  could  optimize  the  distribution  of  various  mesh  parameters  and  lead  to 
superaccurate  methods  which  were  capable  of  delivering  high  accuracy  with  very  few 
degrees  of  freedom  while  capturing  flow  features  that  spanned  many  different  spatial  scales 
of  resolution. 

1.2  Goals  of  the  Project  at  the  Start.  The  research  team  had  some  experience 
developing  hp-type  finite  element  methods  for  simple  heat  conduction  and  viscous  flow 
problems  at  the  start  of  this  project.  With  this  background  and  with  the  possibility  of  high 
dividends  if  the  problem  of  developing  high-level  hp-data  structures  could  be  adequately 
handled,  the  project  was  initiated  and  had  as  its  initial  goals  the  following: 
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1.  Data  structures.  The  development  of  new  finite  element  data  structures  that 
allow  for  arbitrary  mesh  sizes  by  h-refinement  and  arbitrary  spectral  orders  by  p- 
enrichment.  The  data  structure  would  have  to  be  sufficiently  general  to  be  used  in 
modeling  not  only  the  fluid  but  also  the  structure  and  the  interface,  and  would  be  able  to 
resolve  multiple-scale  features  of  flow  phenomena  through  the  use  of  hierarchical  spectral- 
type  approximations. 

2.  High-order  temporal  schemes.  Since  the  time-dependent  case  was  to  be  also 
modeled,  it  was  necessary  to  develop  very  high-order  temporal  approximation  schemes  so 
that  the  order  of  the  approximation  in  time  would  balance  the  corresponding  higher-order 
approximation  in  space. 

3.  Elastic  structure  modeling.  The  same  data  structure  used  for  modeling  the  fluid 
and  the  fluid  interface  was  to  be  configured  so  as  to  be  applicable  to  large-scale  structural 
analysis  as  well. 

4.  New  high-order  methods  for  incompressible  viscous  flow.  At  the  start  of  this 
project,  no  adaptive  h-p  methods  had  ever  been  developed  for  large-scale  incompressible 
flow  problems.  A  major  goal  was  to  also  develop  adaptive  strategies  that  would 
automatically  produce  accurate  results  at  a  specified  scales. 

5.  Fluid-structure  interface.  Various  techniques  for  modeling  the  motion  of  a  solid 
interface  in  viscous  fluid  had  to  be  re-explored  and  adapted  to  these  new  approaches.  This 
involved  a  new  look  at  such  classical  methods  as  arbitrary  Lagrange-Euler  (ALE)  methods, 
etc. 

6.  Iterative  and  flow  solvers.  It  was  realized  at  the  onset  for  large-scale  problems 
many  of  the  direct  methods  for  solving  large  linear  systems  would  not  be  applicable  to  the 
subject  class  of  problems.  A  parallel  study  on  iterative  solvers  for  large-scale  h-p  systems 
was  also  initiated.  Among  these  solvers  were  those  applicable  to  unsymmetrical  systems 
since  these  occur  in  Navier-Stokes  simulations. 

7.  A  posteriori  error  estimation.  The  optimal  control  process  and  hp-adaptive 
schemes  are  based  on  the  concept  of  manipulating  mesh  parameters  so  as  t  j  control  the 
numerical  error.  Consequently,  very  reliable  techniques  for  a  posteriori  error  estimation  are 
needed  to  drive  these  techniques.  At  the  start  of  this  project,  a  posteriori  error  estimates 
were  available  only  for  relatively  simple  linearly  elliptic  problems  and,  except  for  some 
very  special  cases,  none  had  been  developed  and  tested  for  Navicr-Stokes  equations  or  for 
general  fluid  structure  interaction  problems.  The  development  of  reliable  error  estimators 
schemes  was  identified  as  an  important  goal  in  this  research  effort. 
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2.  Technical  Approach 

In  the  early  stages  of  this  project,  attention  was  focused  on  three  of  the  tasks  listed 
in  the  previous  section: 

•  Development  of  h-p  data  jiructures  and  adaptive  strategies 

•  Development  of  high-order  solvers  for  the  incompressible  Navier-Stokes 
equations 

•  Implementation  of  error  estimation  techniques  that  could  be  applicable 

Our  early  work  on  data  structures  met  with  considerable  success.  A  number  of 
papers  were  published  on  the  subject,  and  the  techniques  have  been  the  subject  of 
presentations  at  a  number  of  national  and  international  meetings,  and  extensions  of  the  idea 
are  even  being  used  in  one  or  two  commercial  software  packages  at  the  present.  This  data 
structure  was,  as  far  as  is  known  by  the  author,  the  only  data  structure  at  the  time  that 
allowed  a  genuine  h-p  adaptive  strategy  in  which  the  mesh  size  h  and  the  local  order  p  are 
actually  treated  as  free  parameters.  The  first  such  data  structure  was  developed  by  the 
author  and  his  Ph.D.  students  and  is  summarized  in  Appendix  A,  which  is  excerpted  from 
a  three-part  paper  on  h-p  data  methods,  written  in  1989-90.  Some  of  this  work  formed  the 
basis  of  the  doctoral  dissertation  of  W.  Rachowicz  and  was  developed  jointly  by  the 
author,  L.  Demkowicz,  and  W.  Rachowicz. 

The  development  high-order  adaptive  flow  solvers  for  the  Navier-Stokes  equations 
presented  a  long  series  of  difficulties.  Not  only  is  it  necessary  to  develop  robust  and 
accurate,  high-order  time  marching  schemes,  but  also  notorious  difficulties  with 
approximating  pressures  in  incompressible  flows  had  to  be  dealt  with.  Here  reference  is 
made  to  the  so-called  LBB  condition  in  which  approximations  of  the  velocities  cannot 
necessarily  be  made  independently  of  the  pressures.  A  poor  pressure  approximation  can 
lead  to  unstable  oscillations  and,  for  higher  order  h-p  schemes,  these  oscillations  could  be 
devastating. 

For  these  reasons,  a  number  of  new  high-order  methods  were  first  investigated. 
The  first  among  these  were  the  methods  of  characteristics.  These  schemes,  based  on  rather 
classical  ideas  of  characteristics,  are  used  successfully  in  the  European  aerospace  industry 
to  solve  subsonic  flow  problems.  They  are  inherently  low-order  schemes,  since  they 
involve  effectively  a  splitting  of  the  convection  and  diffusion  steps.  However,  these 
methods  were  easy  to  implement  and  were  among  the  first  to  be  coded  and  tried  for  two- 
dimensional  cases  employing  our  newly  developed  h-p  data  structure. 

Initial  results  were  deceptively  good.  Quite  good  results  for  transient  two- 
dimensional  incompressible  flow  problems  were  initially  obtained.  A  lengthy  series  of 
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studies  of  techniques  for  pressure  approximation  was  also  initiated.  It  was  finally  decided 
that  a  penalty  type  approximation  for  pressures  would  be  most  appropriate  for  the  h-p  data 
structure.This  involved  attempts  at  developing  a  strategy  that  would  always  guarantee 
stable  pressure  approximation.  As  a  general  rule,  this  cannot  be  done  with  exterior  penalty 
formulations  unless  reduced  integration  is  used.  The  subject  of  reduced  integration  was 
studied  extensively  by  the  author  from  1979  to  1982,  but  not  for  h-p  approximations.  This 
represented  new  ground  and  proved  to  be  an  exceedingly  difficult  issue.  Some  theoretical 
results  were  developed  including  a  stability  theorem  that  showed  this:  if  a  Gaussian 
integration  scheme  of  the  order  r  is  required  to  integrate  the  pressure  approximation  exactly, 
then  the  use  of  the  scheme  of  order  r—2  will  always  result  in  a  stable  pressure  (i.e.,  the 
satisfaction  of  the  LBB  condition).  Mathematical  proof  of  this  condition  for  the  two- 
dimensional  Stokes  problem  was  developed  and  numerous  numerical  experiments  were 
done.  This  work  was  never  published,  because  some  of  the  numerical  experiments  led  to 
data  in  contradiction  to  some  of  the  theoretical  results.  Moreoever,  the  penalty 
approximation  is  notoriously  difficult  to  stabilize.  It  was  finally  decided,  after  a  number  of 
months,  that  alternative  techniques  for  pressure  approximation  should  be  investigated. 

In  addition  to  the  method  of  characteristics,  several  new  high-order  methods  of 
characteristics  were  developed.  These  were  coded,  experiments  were  done  on  a  number  of 
two-dimensional  cases,  and  these  methods  were  finally  discarded  as  being  expensive  and 
often  unreliable. 

In  the  literature  of  the  '80s,  considerable  discussion  was  given  to  Galerkin  least- 
square  methods  and  least-square  methods  for  Navier-Stokes  equations.  A  study  of  these 
techniques  was  also  undertaken.  These  were  also  ultimately  discarded  because  they  were 
expensive,  basically  nonlinear  in  character,  and  often  exhibited  very  high  artificial  viscosity 
that  tended  to  obliterate  important  physical  features  of  the  solution. 

More  recently,  a  detailed  study  of  high-order  pressure-based  schemes  was 
undertaken.  These  schemes  originate  from  a  1968  paper  of  A.  Chorin  and  are  perhaps  the 
most  popular  methods  in  use  today  for  solving  the  incompressible  Navier-Stokes 
equations.  They  are  based  on  the  idea  of  first  introducing  a  momentum  step  in  which  the 
velocity  is  advanced  with  or  without  a  pressure  term  in  the  momentum  equation.  Then  a 
Poisson  equation  for  the  pressure,  or  a  pressure  like  variable,  is  solved  and  used  to  correct 
the  velocities.  These  schemes  are  rich  in  mathematical  history,  having  been  studied  by  not 
only  Chorin  but  also  R.  Temam  and  others.  They  are  also  the  basis  of  a  number  of 
commercial  codes. 

Despite  their  popularity,  pressure-based  schemes  are  essentially  first  order  accurate 
in  time.  There  have  been  a  number  of  papers  that  have  attempted  to  extend  these  to  second 
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order  schemes,  but  they  have  been  successful  only  for  special  types  of  boundary 
conditions.  Investigators  of  the  present  project  spent  considerable  effort  in  trying  to 
develop  high  order  pressure  based  schemes.  These  involve  the  use  of  implicit  Runge-Kutta 
methods  for  the  momentum  step  and  the  use  of  h-p  adaptive  solvers  for  the  pressure  step. 
To  date,  it  has  proved  to  be  impossible  to  get  a  reliable  second-order  scheme.  On  the  other 
hand,  we  did  develop  a  robust  first-order  scheme  that  functions  quite  well  on  adaptive  h-p 
meshes.  While  low-order  accurate  in  time,  the  scheme  nevetheless  formed  the  basis  of  our 
numerical  experimentation  up  until  the  conclusion  of  the  project. 

On  another  front,  modeling  of  the  elastic  structure  and  the  interface  was  also 
investigated.  A  high-order  adaptive  h-p  algorithm  was  developed  for  linear  elasticity 
problems.  This  scheme  was  also  used  in  several  experiments  to  model  plates  and  shells 
with  varying  degrees  of  success.  Also,  some  effort  was  spent  on  exploring  methods  for 
modeling  moving  interfaces.  The  result  of  these  studies  was  a  technique  that  was  basically 
an  Arbitrary  Lagrange-Eulerian  scheme,  but  it  did  include  the  capability  of  handling  h- 
adaptivity.  Plans  were  to  eventually  upgrade  this  interface  scheme  to  h-p  methods  as  well, 
but  this  phase  of  the  study  was  never  completed  because  of  focus  of  the  project  on  other 
issues. 

In  addition  to  the  work  on  data  structures  and  flow  solvers  for  h-p  adaptivity,  a 
great  deal  of  the  project  focused  on  the  development  new  a-posteriori  error  estimators.  The 
ability  to  quickly  estimate  the  local,  element-wise  approximation  error  during  the  evolution 
of  a  calculation  is  a  key  feature  of  any  successful  adaptive  scheme.  At  the  initiation  of  this 
project,  the  available  technology  for  error  estimation  was  rather  primitive  being  primarily 
confined  to  ad  hoc  methods  for  h  adaptivity  based  on  estimates  of  gradients  of  density  or 
pressure;  more  rigorous  techniques  were  successfully  applied  only  to  very  limited  classes 
of  one-dimensional  elliptic  problems.  During  early  stages  of  the  project,  considerable 
progress  was  made  on  the  development  of  the  error  residual  method,  a  concept  first 
introduced  by  the  present  team  in  the  early  '80s.  This  method  was  further  developed  and 
applied  to  error  estimation  in  linear  elasticity  problems  and  finally  to  Navier-Stokes 
equations.  A  paper  detailing  some  of  the  earlier  work  on  a  posteriori  error  estimation  for 
linear  elliptic  problems  was  produced  early  in  the  project  and  eventually  implemented  and 
applied  to  the  Navier-Stokes  equations. 

More  recently,  a  new  variant  of  the  error  residual  method  was  derived  which  is 
based  on  so-called  flux  balancing.  It  was  discovered  during  the  course  of  the  work  that  the 
quality  of  an  error  estimator  depends  very  much  on  the  ability  to  approximate  error  fluxes  at 
the  boundaries  of  finite  elements.  Prior  to  this  work,  no  error  estimation  technique  was 
known  that  would  function  well  for  arbitrary  h-p  meshes.  The  basic  idea  behind  these 
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techniques  is  this:  the  error  in  a  gridcell  depends  upon  the  local  (element-wise)  residual 
which  is  a  calculable  function  representing  the  degree  to  which  the  exact  solution  fails  to 
satisfy  the  governing  equations.  In  addition  to  this,  the  error  is  influenced  by  the  boundary 
conditions  on  an  element,  the  so-called  flux  of  error  from  one  cell  into  another.  The  idea 
behind  the  error  residual  method  is  to  compute  element  residual  and  set  up  a  local  boundary 
value  problem  which  is  solved  quickly  for  a  running  account  of  the  local  error.  If  the 
original  problem  is  solved  on  a  mesh  for  which  polynomial  approximation  of  degree  p  are 
used  in  a  gridcell,  the  error  estimate  was  generally  done  by  estimating  the  error  with 
polynomials  of  degree  p  or  higher.  Generally,  a  higher  p  was  required  for  the  error 
estimator  in  regions  of  singularities,  such  as  in  separation  points,  changes  in  boundary 
conditions  or  at  point-line-singularities.  It  was  discovered  that  for  odd-order 
approximations,  a  key  factor  in  the  quality  of  the  estimate  was  the  balance  of  the  error 
fluxes  and  a  new  technique  was  developed  to  resolve  this  issue.  This  work,  recently 
accepted  for  publication  in  Numerische  Mathematik,  provides  a  significant  advance  to  the  a 
state-of-the-art  posteriori  error  estimation  and  has  proved  to  be  successful  in  certain  two- 
dimensional  calculations  of  unsymmetrical  elliptic  systems. 

Finally,  some  additional  comments  should  be  made  concerning  iterative  solvers  and 
direct  solvers.  A  portion  of  the  effort  was  spent  investigating  various  types  of  iterative 
solvers  with  the  understanding  that  large-scale  fluid  structure  interaction  problems  cannot 
be  efficiently  handled  by  standard  direct  solvers  or,  induced,  by  direct  solvers  in  general. 
This  part  of  the  work  did  not  lead  to  any  new  methodologies  but  rather  an  assimilation, 
comparison,  evaluation,  and  implementation  of  a  number  of  existing  schemes.  Finally,  in 
much  of  the  work  we  continued  to  use  a  direct  frontal  solver  because  of  its  robustness  and 
because  the  purpose  of  the  early  numerical  experimentation  was  not  necessarily  to  solve 
specific  large-scale  problems,  but  to  explore  the  efficiency  and  functuality  of  new  methods. 
However,  we  did  implement  and  employ  various  iterative  schemes  such  as  the  GMRES 
solver,  Jacobi-Block  iteration  and  other  techniques.  This  work  should  prove  to  be  valuable 
later  when  these  techniques  are  extended  to  large-scale  problems  and  to  multi-processor 
environments. 

3.  Progress 

3.1  Successes  and  Failures.  In  summarizing  the  approaches  outlined  in  the  previous 
section,  one  can  point  to  a  number  of  successes  and  a  number  of  failures  due  to  "blind 
alleys"  pursued  during  the  project.  The  major  successes  center  around  three  areas:  the 
development  of  new  h-p  adaptive  strategies  and  data  structures,  the  development  of  new 
methods  for  solving  the  incompressible  Navier-Stokes  equations  and  linear  elasticity 
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equations  using  h-p  adaptive  solvers,  and  the  development  of  a  new  family  of  a  posteriori 
error  estimators. 

As  noted  above,  a  number  of  the  initial  approaches  led  to  "  blind  alleys."  These 
included  in  particular  the  lengthy  studies  of  penalty  approaches  to  h-p  adaptive  schemes  for 
Navier-Stokes,  the  development  of  high  order  methods  of  characteristics,  which  while 
moderately  successful  did  not  appear  to  provide  a  robust  foundation  for  further  work,  the 
exploration  of  Galerkin  and  least-squares  based  methods  for  primitive  variable 
approximations  of  incompressible  flow,  and  generally  the  development  of  high-order 
temporal  schemes  for  Navier-Stokes  equations.  Indeed,  this  latter  area  remains  open  and 
represents  one  critical  area  in  which  considerable  additional  work  needs  to  be  done.  At 
present,  our  experimental  codes  function  on  two  classes  of  algortihms:  1)  pressure  based 
schemes,  which  employ  a  singly-implicit  Runge-Kutta  method  for  the  momentum  step,  and 
the  full  high-order  two-stage  Runge-Kutta  methods  for  the  entire  Navier-Stokes  equations. 
There  is  no  doubt  that  these  schemes  can  be  significantly  improved  with  additional  study. 

A  summary  paper,  including  much  of  our  work  on  data  structures  and  some  of  our 
experimentation  with  high-order  schemes  in  attached  in  Appendix  B  of  this  repon. 

4.  Changes:  Final  Approaches  and  Final  Goals 

The  success  with  the  h-p  adaptive  schemes  throughout  the  project  made  it  clear  that 
this  is  indeed  a  viable  approach  for  large-scale  fluid-structure  interaction  calculations.  The 
data  structures  originally  developed  in  this  project  have  undergone  a  number  of 
evolutionary  changes  and  enough  information  has  been  gathered  to  produce  new  versions 
of  these  schemes  which  should  be  highly  effective.  The  error  estimation  techniques 
developed  seem  to  be  quite  promising,  but  additional  work  is  needed  to  define  their 
limitations  and  to  extend  them  to  truly  coupled  problems  of  fluid  structure  interaction  and  to 
rigorous  methods  of  the  full  Navier-Stokes  equations.  The  final  stages  of  the  project 
focused  primarily  on  the  Navier-Stokes  simulation  since  the  methods  developed  seem  to 
work  quite  well  for  modeling  structural  response. 

5.  Relationship  to  Other  Work 

It  should  be  noted  that  a  number  of  the  key  components  of  the  present  project 
overlapped  with  on-going  work  supported  in  a  parallel  ONR  project  on  structural  acoustics. 
Indeed,  much  of  the  experiments  developed  on  error  estimation  and  on  h-p  data  structures 
proved  to  be  invaluable  in  developing  similar  procedures  for  BEM  techniques  for  structural 
acoustics  problems.  In  those  problems,  the  primary  modeling  technique,  however,  was 
based  on  boundary  element  methods  and  somewhat  different  data  structures  and  error 


estimation  techniques  applied.  On  the  other  hand,  experiments  developed  in  both  the  fluid 
structure  interaction  effort  and  in  the  structural  acoustics  effort  proved  to  be  beneficial  in 
each  of  the  projects. 

6.  Problems  Perceived 

In  future  work,  the  direction  toward  more  useful  results  is  quite  clear.  Additional 
work  on  developing  robust  high-order  solvers  that  function  on  high-order  h-p  meshes 
needs  to  be  done.  Continued  work  on  pressure  based  schemes  is  intriguing  and  should  be 
continued. 

During  the  first  year  of  this  project  we  developed  what  appears  to  be  the  first  h-p 
adaptive  strategy  and  has  been  used  successfully  in  the  literature.  This  refers  to  a  technique 
applicable  to  arbitrary  h-p  meshes,  which  when  coupled  with  estimates  of  the  local  error, 
allows  for  a  decision  to  be  made  whether  to  refine  or  coarsen  the  mesh,  or  enrich  or  de- 
enrich  the  spectral  order.  This  is  also  a  key  to  the  success  of  h-p  adaptive  methods.  Much 
additional  work  needs  to  be  done  on  this  subject.  In  ’•ecent  months,  a  number  of  new 
adaptive  strategies  have  been  developed  and  have  been  implemented.  These  show  in  some 
cases,  exponential  convergence  of  the  error  versus  CPU  time,  as  opposed  to  degrees  of 
freedom,  and  are  regarded  as  significant  advances  in  this  subject.  These  methods  also  lend 
themselves  to  parallelization,  and  the  implementation  of  these  methods  on  parallel 
processers  is  suggested  as  a  fundamental  and  important  component  of  any  additional  work 
that  is  done  on  this  subject. 

It  is  believed  that  the  basic  components  of  an  effective  adaptive  strategy  for  quite 
large  problems  in  fluid  structure  interaction  is  in  hand.  Future  efforts  should  include 
significant  work  on  parallelization  and,  the  development  of  load-balancing  techniques  for 
adaptive  strategies.  Indeed,  while  h-p  adaptive  methods  or  serial  computers  could  very 
well  represent  the  reduction  of  computational  times  by  one  to  two  orders  of  magnitude,  the 
implementation  of  these  schemes  on  multi-processor  platforms  could  very  well  provide  an 
additional  order  of  magnitude  or  more  in  increases  in  efficiency. 

As  to  major  problems  that  are  perceived,  a  significant  problem  may  have  to  do  with 
the  stability  and  conditioning  of  h-p  methods.  With  the  increase  in  spectral  order,  one  also 
finds  an  increase  in  the  condition  number  of  the  underlying  matrices.  Progress  has  been 
made  to  develop  pre-conditioning  methods  to  control  the  stability  of  these  schemes,  but  this 
is  still  a  young  area  of  research  and  much  additional  work  needs  to  be  done.  The 
parallelization  and  preconditioning  of  adaptive  h-p  methods  would  clearly  represent  a 
fundamentally  important  area  of  study  that  needs  to  be  done  prior  to  the  extension  of  these 
technologies  to  significant  practical  problems  in  modeling  fluid- structure  interactions. 


7.  Documentation 


The  research  outlined  above  led  to  the  publication  of  a  number  of  articles,  reports, 
and  papers.  Also  several  oral  presentations  of  the  work  were  made  at  national  and 
international  conferences  and  symposia.  A  summary  of  these  documents  and  lectures  is 
given  below. 
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APPENDIX  A 


Constrained  Approximation:  A  General  Approach  to 
h-p  Adaptive  Data  Structures 


In  the  body  of  this  report,  refrence  is  made  to  various  hp-data  structures  developed  and 
studied  during  the  course  of  this  work.  This  appendix  is  provided  to  furnish  some  detail 
on  one  such  data  structure  for  completeness  and  for  reference  purposes.  The  particular 
formulation  outlined  here  follows  a  paper  by  Demkowicz,  Oden,  Rachowicz,  and  Hardy 
[CMAME,  Vol.  77,  pp.  79-112]  and,  while  not  the  most  sophisticated  version  in  use, 
provides  all  of  the  key  details  and  features  needed  to  understand  the  basic  strategy. 


1  h-  and  p- Adaptivity.  Regular  and  Irregular  Meshes 


Let  fi  be  an  open  bounded  domain  in  SP',  n  =  2, 3  with  a  sufficiently  regular  boundary  3f2. 
In  what  follows,  we  shall  restrict  ourselves  to  a  class  of  problems  that  can  be  formulated  in 
the  following  abstract  form: 


Find  u  ^  X  such  that 


B{u,v)  =  Liv)  Vv  e  X  j 


(1.1) 


Here: 


X  =  X  X  X  ■  • '  X  X{Tn  times)  (1-2) 

where  X  a  subspace  of  77^(17),  the  Sobolev  space  of  first  order,  is  a  bilinear  form  on 

X  X  AC  of  the  following  form 


B{u,v)  =  ^  Bij(ui,Vj) 

».j=i 


(1.3) 


where  5,j(-,  •)  are  bilinear  forms  of  scalar-valued  arguments  and  L{-)  is  a  linear  form  on  X 
of  the  form 

L(v)  =  ■£  Liivi)  (1.4) 
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Finite  Element  Approximation 


We  assume  that  the  domain  can  be  represented  as  a  union  of  finite  elements  A'e,  e  = 
1,. .  ■  ,M.  More  precisely 

_  M _ 

n=[jK.  (1.5) 

e=l 

and 

int  A'e  n  int  Kj  =  (f)  for  e  7^  /  (1-6) 

Each  of  the  elements  K  has  a  corresponding  finite  dimensional  space  of  shape  functions, 
denoted  Xh{K),  for  instance  the  space  of  polynomials  of  order  p.  The  global  finite  element 
space  Xh  consists  of  functions  which,  when  restricted  to  element  K,  belong  to  the  local 
space  of  shape  functions  Xh{K)-  Thus  the  global  approximation  is  constructed  by  patching 
together  the  local  shape  functions  in  the  usual  way. 

We  shall  adopt  the  fundamental  requirement  that  the  global  approximation  must  be  con¬ 
tinuous.  As  we  will  see,  this  requirement  leads  to  the  notion  of  constrained  approximations. 
Formally,  the  continuity  assumption  guarantees  that  the  finite  element  space  X^  is  a  sub¬ 
space  of  and,  with  some  additional  assumptions  if  necessary,  also  a  subspace  of  X. 

The  approximate  problem  is  easily  obtained  from  (2.1)  by  substituting  for  u  and  v  their 
approximations  Uh  and  Vh: 


Find  Uh  S  Xh  such  that  1 


Bhiuh^Vh)  =  Lh{vh)  ^Vh^Xh  j 


(1.7) 


Here 

Xh  =  Xh  X  ...  X  Xh{m  times) 


(1.8) 


which  indicates  that  the  same  approximation  has  been  applied  to  every  component  of  u. 
Bh{-,  •)  and  Lh{-)  denote  approximations  to  the  original  bilinear  and  linear  forms  resulting 
from  numerical  integration. 


Adaptivity 

A  flow  chart  of  a  typical  Adaptive  Finite  Element  Method  (AFEM)  is  shown  in  Fig.  1. 
The  method  consists  of  first  generating  an  initial  mesh  and  solving  for  the  corresponding 
FEM  approximate  solution.  Next,  the  error  is  estimated  in  some  way  and  based  on  this 
(usually  crude)  approximation,  one  adapts  the  mesh,  i.e.,  adds  new  degrees  of  freedom.  The 
approximate  problem  for  the  new  mesh  is  solved  again  and  the  whole  procedure  continues 
until  certain  error  tolerances  are  met.  Obviously,  such  a  procedure  requires  an  estimate  of 
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the  error  over  each  element  and  a  strategy  to  reduce  the  error  by  proper  changes  in  the  mesh 
parameters,  h  and  p. 

In  the  FEM  the  new  degrees  of  freedom  can  be  added  in  two  different  ways:  elements 
may  be  locally  refined  or  their  spaces  of  shape  functions  may  be  enriched  by  incorporating 
new  shape  functions.  As  noted  earlier,  in  the  caise  of  polynomials,  this  is  done  by  increasing 
locally  the  degree  of  polynomials  used  to  construct  the  shape  functions,  the  first  case  being 
an  A-refinement,  and  the  second  case  a  p-refinement.  A  combination  of  both  is  an  (adaptive) 
h-p  FEM.  We  remark  that  the  process  of  increasing  the  local  polynomial  degrees  for  a  fixed 
mesh  size  is  mathematically  akin  to  increasing  the  spectral  order  of  the  approximation  and 
that,  therefore,  we  also  refer  to  h-p  methods  as  “adaptive  spectral-element”  or  “h-spectral'’ 
methods. 

Regular  and  Irregular  Meshes 

As  the  result  of  local  h-refinements,  irregular  meshes  are  introduced.  Recall  (see  [6])  that  a 
node  is  called  regular  ii  it  constitutes  a  vertex  for  each  of  the  neighboring  elements;  otherwise 
it  is  irregular.  If  all  nodes  in  a  mesh  are  regular,  then  the  mesh  itself  is  said  to  be  regular. 
In  the  context  of  two-dimensional  meshes,  the  maximum  number  of  irregular  nodes  on  an 
element  side  is  referred  to  as  the  index  of  irregularity.  Meshes  with  an  index  of  irregularity 
equal  one  are  called  1-irregular  meshes.  The  notion  can  be  easily  generalized  to  the  three- 
dimensional  case.  (See  [7]  and  literature  cited  therein  for  additional  references.) 

In  the  present  work,  we  accept  only  l-irregular  meshes.  In  the  two-dimensional  context, 
this  translates  into  the  requirement  that  a  “large”  neighbor  of  an  element  may  have  no  more 
than  two  “small”  neighbors  on  a  side;  in  the  three-dimensional  case,  the  number  of  neighbors 
sharing  a  side  may  go  up  to  four,  while  the  number  of  neighbors  sharing  an  edge  can  be  no 
more  than  two.  This  is  frequently  called  the  “two-to-one”  rule  (cf.  [7]).  Examples  of  regular 
and  irregular  meshes  are  shown  in  Fig.  2.  There  are  several  practical  and  theoretical  reasons 
to  accept  only  l-irregular  meshes,  especially  in  the  context  of  h-p  methods.  For  a  detailed 
argument,  we  refer  to  [8]. 

Our  restriction  to  l-irregular  meshes  imposes  a  simple  restriction  on  the  way  any  h- 
refinement  can  proceed:  before  an  element  is  refined,  a  check  for  “larger”  neighbors  must  be 
made.  If  any  such  neighbors  exist,  they  must  be  refined  first  and  only  then  can  the  element 
in  question  be  refined. 
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Figure  1:  Typical  flow  chart  of  an  adaptive  method 
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Figure  2:  Exai 
1- irregular  mesl 


Continuity  and  Constrained  Approximation 

The  presence  of  irregular  nodes  makes  the  handling  of  the  continuity  assumption  non¬ 
standard  and  leads  to  the  notion  of  a  constrained  approximation.  As  an  example,  consider  a 
mesh  of  three  rectangular  elements  (quadrilateral  elements  with  biquadratic  shape  func¬ 
tions)  with  standard  Lagrange  degrees  of  freedom,  as  shown  in  Fig.  3.  Clearly,  a  function 
u^^  defined  over  these  elements  need  not  to  be  continuous  across  element  interfaces  due  to 
the  presence  of  two  irregular  nodes,  A  and  B.  For  instance,  the  values  of  on  side  CDE  of 
element  K  are  determined  uniquely  by  (say)  its  values  at  points  C,  D  and  E,  while  those  on 
side  CD  of  element  K2  are  defined  by  specifying  degrees  of  freedom  in  K.  In  order  to  make 
u/,  continuous,  the  value  of  Uh  at  A  from  the  side  of  K2  must  be  forced  to  be  equal  to  the 
value  of  Ufi  at  A  from  the  side  of  Kx ,  which  is  equivalent  to  the  elimination  of  the  degree  of 
freedom  associated  with  point  A  by  enforcing  the  constraint 


Uh{A)  =  auh{C) -f- Puk{D)  + -fUh{E)  (1.9) 


with  proper  coefficients  a,  ^  and  7.  The  situation  is  much  more  complicated  in  the  case  of 
different  orders  of  approximation  within  each  element,  and  we  discuss  it  in  detail  in  Section 
3. 

2  Constrained  Approximation 

In  this  section,  we  develop  the  general  concept  of  constrained  approximations.  We  shall  see 
in  particular  the  impact  of  the  constraints  on  such  basic  ingredients  to  the  FEM  as  element 
stiffness  matrix  and  load  vector  calculations. 

Approximation  on  the  Element  Level 

Let  K  be  a  finite  element  with  the  corresponding  space  of  shape  functions 

The  set  Qfc  of  element  degrees  of  freedom  W/c,  as  usual,  is  viewed  as  a  set  of  linear  functionals 
defined  on  Xh{K).  We  assume  that  the  set  Qk  of  degrees  of  freedom 


Wi,K  •  Xh{K)  — I  i  e  Nk  C  N} 


(2.1) 
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Figure  3;  Example  of  an  unconstrained,  discontinuous  approximation. 
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is  X/i(/<')-unisolvent  (see  [9,  10]),  i.e., 

1°  €  Nfc,  are  linearly  independent 

2°  6  Nk  is  a  dual  space  to  X^iK) 

As  indicated,  we  assign  to  each  linear  functional  G  Qk  an  integer  label  i,  p  — >  and 

denote  by  Nk  the  set  of  such  labels  for  element  K.  The  element  shape  functions  Xj,K  are 
defined  as  a  dual  basis  to  he.. 


<  ^i,K-,Xj,K  >=  ^ij  hj  6  Wa' 


(2.2) 


and  the  finite  element  approximation  Uh  within  the  element  K  is  of  the  form 

=  Zi  <Xi,K 

i€Nk 

where  =  <fijc{uh). 

In  what  follows,  we  restrict  ourselves  to  Lagrange-  and  Hermite-types  of  degrees  of 
freedom.  In  other  words,  we  assume  that  each  of  the  functionals  pi  is  of  the  form 


p  :  u 


^  —  0,1,... 


(2.3) 


where  DxU  denotes  the  k-th  order  differential  of  u  evaluated  at  point  x  (as  usual,  Dj-u  = 
u(®)).  Vectors  at  this  point  denote  arbitrary  vectors  in  n  =  2,3.  Thus,  every 

degree  of  freedom  can  be  identified  with  a  point  x  and  k  vectors  (Later  on,  in 

the  context  of  subparametric  finite  elements,  we  shall  accept  functionals  which  are  linear 
combinations  of  functionals  of  the  form  (3.3).) 

Construction  of  the  Unconstrained  Finite  Element  Space  X* 

We  introduce  the  following  formal  definition  of  the  unconstrained  finite  element  space  Xh. 
A  function  u/,:  fl  — belongs  to  Xh  if  and  only  if  the  following  two  conditions  are  satisfied: 

1.  The  restriction  of  Uh  to  an  element  K  belongs  to  the  local  space  of  shape  functions, 
i.e., 

uhIk  G  Xh{K)  VA'  (2.4) 
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2.  For  every  pair  of  elements  Ke  and  K /  and  corresponding  pair  of  degrees  of  freedom 

such  that  and  tpiCf  are  defined  by  (3.3)  through  the  same,  common  point  x  and 
vectors  degrees  of  freedom  take  on  the  same  value  on  u/,  when 

restricted  to  /fg  or  Kj  respectively,  i.e., 


9KA^h\Kt)  =  'pK,Wh\Ki) 


(2.5) 


Note  that  the  elements  of  Xh.  need  not  be  continuous  (see.  Fig.  3). 

Global  Degrees  of  Freedom 

Due  to  the  const  iction  of  the  space  Xh,  we  can  introduce  the  global  degrees  of  freedom 
identified  by  points  (nodes)  x  and  vectors  We  define  the  linear  functional  $  on 

Xh, 

^:Xh'^R,H^h)^<PK{nh\K),  (2.6) 

where  K  is  an  element  with  the  corresponding  degree  of  freedom  identified  with  point  x  and 
vectors  . . .  ,^h-  Note  that  due  to  the  definition  of  Xh,  the  global  degrees  of  freedom  are 
well  defined. 

The  Unconstrained  Base  Functions 

The  unconstrained  base  functions  e,-  are  introduced  as  a  dual  basis  to  the  space  of  the  global 
degrees  of  freedom,  i.e., 

(2.7) 

Note  that  e,-  may  be  discontinuous. 

Construction  of  the  Constrained  Finite  Element  Space  Xh 

At  this  point,  somewhat  arbitrarily,  we  divide  all  global  degrees  of  freedom  into  two  subsets: 
active  and  constrained.  We  use  the  notation. 
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•  the  set  of  indices  aissigned  to  active  degrees  of  freedom  and 

•  the  set  of  indices  assigned  to  constrained  degrees  of  freedom 

By  an  “active”  degree  of  freedom,  we  mean  one  of  a  set  of  linearly  independent  functionals 
that  define  the  parameters  associated  with  a  global  stiffness  matrix  for  the  problem  at  hand; 
“constrained”  degrees  of  freedom  are  linear  combinations  of  active  degeres  of  freedom  defined 
by  constraints  associated  with  element  connectivity. 

We  assume  that  for  each  constrained  degree  of  freedom  i  €  there  exists  a  set  I{i) 
of  corresponding  active  degrees  of  freedom,  I{i)  C  iV“,  and  a  vector  j  G  I{i),  such  that 
the  following  equality  holds: 

^i{uh)  =  (2.8) 

i6/(.) 

We  now  introduce  the  constrained  finite  element  space  Xh.  as 


X,.  =  u,,  6  X,,  I  $.(u0  =  E  Vf  e  N‘=  \  (2.9) 

I  j€/(.)  J 


Assuming  that  the  constraints  are  linearly  independent,  we  see  that  Xh  is  dual  to  the 
space  spanned  by  only  active  degrees  of  freedom.  As  usual,  we  define  the  base  functions  ej, 
j  €  W“,  as  a  dual  basis  to  the  set  of  active  degrees  of  freedom: 


ej  e  Xh,  <  >=  Sij  i,j  G 


(2.10) 


Though,  at  this  point,  the  choice  of  constrained  degrees  of  freedom  is  arbitrary,  we 
implicitly  assume  that  with  the  proper  choice  of  constraints  the  resulting  finite  element 
space  Xh  consists  of  only  continuous  functions. 
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Relation  Between  Unconstrained  and  Constrained  Base  Functions 

Let  Uh  be  an  arbitrary  function  belonging  to  X^-  Then  Uh  must  be  of  the  following  form: 

Uh  = 

(2.11) 

=  ^  w.e,  +  RjkUkej 

ieN‘^  jeN-^  k€l{j) 

Introducing  for  every  i  €  iV“  the  set 

S(i)  =  {j  6  iV‘=  1  i  e  I(J)}  (2.12) 

we  rewrite  in  the  form 

Uh  =  Y  m  UhRjkej 

i^N<^  keN‘‘  j€S{k) 

=  n  (e.  +  Z) 

ieN-  \  jes{i)  ) 

We  claim  that  functions 

Ci  =  e,-  T  ^  "  RjiCj  i  G  W** 

ie5(i) 

form  the  dual  basis  to  functionals  i  G  iV“.  Indeed 

<  $j,ei  >=<  >  +  Y  ^ki<  >=  6,j  (2.15) 

ife€S(i) 

since  S{i)  C  N^. 

Calculation  of  the  Global  Load  Vector  and  Stiffness  Matrix 

For  simplicity,  we  restrict  ourselves  to  the  case  of  a  single  equation.  In  the  case  of  systems, 
the  same  procedure  is  applied  for  every  linear  form  (2.7)  and  bilinear  form  (2.4). 


(2.13) 

(2.14) 
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Substituting  (3.14)  into  both  sides  of  (2.10),  we  obtain  formulas  for  the  load  vector  and 
stiffness  matrix: 

Lh{^i)  =  Lh{ei) ^  RkiLhick)  (2-16) 

kes(i) 


+ 

keS(i) 

+  ^  RiiBh{ei,ei) 

tes(j) 


(2.17) 


+  B.kiRijBh{ek.,et) 

kes(i)  les(j) 


Element  Level  Revisited — Modified  Element  Stiffness  Matrix  and  Load  Vec¬ 
tor 

Consider  an  element  K.  Let  iV“(/<')  and  N'^{K)  denote  sets  of  indices  corresponding  to 
active  and  constrained  degrees  of  freedom  for  the  element.  Assuming  that,  as  usual,  the  load 
vector  and  stiffness  matrix  are  calculated  by  summing  up  the  contributions  of  all  elements, 
i.e., 

Lhjuh)  =  y^^Lh.K(uh\K)  (2.18) 

K 

and 

Bh{UhiVh)  =  y~'.Bh.Kiuh\K,Vh\K)  (2.19) 

K 

we  arrive  at  the  practical  issue  of  how  to  calculate  the  contributions  of  element  K  to  the 
global  load  vector  and  stiffness  matrix.  Toward  resolving  this  question,  we  introduce 

•  the  usual  element  load  vector 


Kk  =  LH,K{Xi,K)  i  e  N^{K)  U  N‘^{K)  (2.20) 
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•  the  element  stiffness  matrix 


Bi:,K  =  BHAXi,K.Xi,K)  ij  e  N^{K)  U  N%K)  (2.21) 


•  the  set  of  associated  active  degrees  of  freedom  for  element  K 


A^(/0  =  iV“(/0  u  (J  I{j)  (2.22) 

jeNc(K) 


Notice  that  the  two  sets  on  the  right-hand  side  of  (3.22)  need  not  be  disjoint! 

•  the  element  contribution  to  the  global  load  vector  (modified  element  load  vector) 


bi,K  =  W(eik)  i  6  NiK)  (2.23) 


•  the  element  contribution  to  the  global  stiffness  matrix  (modified  element  stiffness  ma¬ 
trix) 

Bij,K  =  Bfi^f({ei\K,e.j\i()  i,j  G  N{K)  (2.24) 
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Algorithms  for  the  Calculation  of  bij^  and 


The  load  vector  and  stiffness  matrix  for  the  constrained  degrees  of  freedom  are  now  computed 
by  transformations  embodied  in  the  following  sample  algorithms: 

FOR  i  e  N{K) 

hi,K  =  0 

ENDFOR 

FOR  i  e  N\K)  U  N^(K) 


ENDFOR 


IF  i  e  iV“(A^)  THEN 
bi,K  =  bi^K  +  bi^K 
ENDIF 

IF  i  e  N-^iK)  THEN 

FOR  k  e  I{i) 

[bk,K  =  bk,K  +  tiikbi,K] 
ENDFOR 

ENDIF 
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The  Algorithm 


FOR  i,j  €  N{K) 

Bij,K  =  0 

ENDFOR 

FOR  ij  e  N'^iK)  U  N^{K) 

IF  i  e  N‘^{K)  AND  j  e  THEN 

Bij,K  =  Bij^K  +  Bij^K 

ENDIF 

IF  i  e  N^K)  AND  j  e  N‘^{K)  THEN 

FOR  k  6  I(i) 

Bkj,K  =  Bkj,K  +  RikBkj,K 

ENDFOR 

ENDIF 

IF  i  e  Ar“(A')  AND  j  G  NV<)  THEN 

FOR  /  G  I{j) 

Bii,K  =  Bii^K  +  B.jiBii,K 

ENDFOR 

ENDIF 

IF  i  G  N\K)  AND  j  G  N‘{K)  THEN 

FOR*G/(0,^e/(i)  _ 

Bki,K  =  Bki,K  +  RikRjiBki,K 

ENDFOR 

ENDIF 

ENDFOR 

3  An  h-p  Adaptive  Finite  Element  Method 

This  section  is  devoted  to  a  presentation  of  an  h-p  Adaptive  Finite  Element  Method  which 
provides  for  instantaneous  h-  and  p-refinements  and  unrefinements.  For  simplicity,  the  dis¬ 
cussion  is  restricted  to  the  two-dimensional  case. 

We  shall  adopt  the  following  assumptions: 

•  the  initial  mesh  is  (topologically)  a  portion  of  a  regular,  rectangular  grid  in 

•  only  1-irregular  meshes  are  accepted  for  all  h-refinements; 


•  the  local  order  of  approximations  may  differ  in  each  element; 

•  the  approximation  must  be  continuous. 

One-Dimensional  Hierarchical  Element  of  an  Arbitrary  Order  p 

We  begin  our  discussion  with  the  definition  of  the  standard  one-dimensional  hierarchical 
master  element  defined  on  the  interval  [—1,  1]. 

The  element  degrees  of  freedom  are  identified  with  three  nodes:  two  endpoints  and  the 
midpoint  x  =  0.  If  the  element  is  of  the  first  order,  only  two  degrees  of  freedom  are  present 
—  these  are  the  nodal  values  at  x  =  — 1  and  x  =  1.  Starting  with  p  =  2,  new  degrees  of 
freedom  are  added  which  are  of  the  form 


V?p(«) 


1  d^u 

\p  dxP 


(0) 


(3.25) 


where  Ap  is  a  scaling  factor. 

The  corresponding  shape  functions  are  illustrated  in  Fig.  4. 

Definition  of  the  Hierarchical  Square  Master  Element 

Setting  K  =  [—1,1]  x  [—1,1],  we  define  the  space  of  shape  functions  as 


Xt,{K)  =  Q^{K) 


(3.26) 


where  denotes  the  space  of  polynomials  up  to  p-th  order  with  respect  to  each  of  the 
variables  separately. 

The  degrees  of  freedom  are  defined  cis  follows: 

•  function  values  at  four  vertices: 


w(-l»-l)»«(l,-l),w(l,l),u(-l,l) 


(3.27) 
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Figure  4:  One-dimensional  hierarchical  master  element.  Degrees  of  freedom  and  correspond¬ 
ing  shape  functions. 
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•  tangential  derivatives  (up  to  a  multiplicative  constant)  up  to  p-th  order  associated  with 
the  midpoints  of  the  four  edges: 


d^iL 

^.-^(1,0)  *  =  2,...,p 

A:  =  2,...,p 
‘  =  ' . " 


mixed  order  derivatives  associated  with  the  central  node 


(3.28) 


(3.29) 


One  can  easily  see  that  the  space  of  the  shape  functions  Q’’{K)  is  a  tensor  product  of 
PP[— 1, 1]  with  itself  and  the  degrees  of  freedom  just  introduced  are  simply  the  tensor  products 
of  the  degrees  of  freedom  for  the  one-dimensional  element.  More  precisely,  if  u  6  Q’’,  then 
u  is  of  the  form 

k 

where  Vk,  Wk  €  P'’(— 1, 1),  and  each  of  the  degrees  of  freedom  can  be  represented  in  the  form 

=  '^u'’{(piiS)^j){vk<S)Wk)  (3.31) 

k 

This  in  particular  implies  that  the  corresponding  shape  functions  can  be  identified  with  the 
tensor  products  of  one-dimensional  shape  functions,  which  are  of  the  form 


Xii^)Xj{y)  i,j  =  0,l,...,p 


(3.32) 


For  i,j  =  0, 1,  we  have  the  usual  bilinear  element  with  four  nodal  degrees  of  freedom. 
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Master  Element  of  an  Enriched  Order 

Let  K  be  an  element  of  p-th  order.  By  adding  additional  shape  functions  corresponding 
to  the  (p  +  l)th  order,  together  with  the  corresponding  degrees  of  freedom,  we  obtain  a 
well-defined  finite  element  whose  space  of  shape  functions  includes  Q^.  In  particular,  by 
adding  all  the  additional  degrees  of  freedom  corresponding  to  the  (p  +  l)th  order,  we  pass 
to  the  element  without  modifying  the  existing  shape  functions  and  degrees  of  freedom 
of  QP.  Equivalently,  we  can  begin  an  adaptive  process  with  element  and  eliminate  some 
degrees  of  freedom  to  pass  to  an  incomplete  element  of  a  lower  order. 

Subparametric  Hierarchical  Elements 

Consider  the  master  element  of  (possibly)  an  incomplete  order  p.  Even  though  p  can  be 
arbitrarily  large,  the  element  may  be  only  Q^-complete,  which  means  that  some  of  the  nodes 
may  be  missing.  An  example  of  such  an  element  is  presented  in  Fig.  5.  An  arbitrary 
location  of  the  seven  nodes  in  the  plane  (x,y)  determine  uniquely  a  map  T  from  the  master 
element  into.ffi^,  the  components  of  T  belonging  to  the  incomplete  space.  More  precisely, 
if  =  1,. . .  ,9  are  the  regular  shape  functions  for  the  nine  node  biquadratic  element,  then 


9 


T{x,y)  = 

iai 


(3.33) 


with  the  assumption  that  og  =  ^(02  +  03)  and  07  =  ^(03  +  04). 
We  have  the  classical  definition  of  the  subparametric  element. 


K  =  T{K) 


(3.34) 


with  the  space  of  test  functions  defined  as 


Xh{K)  =  {u  =  u-t-'  \  uex,{K)} 


(3.35) 
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Figure  5:  Concept  of  the  subparametric  element. 
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and  the  degrees  of  freedom 


<  <p,u  >=<  (p,u  >  where  u  =  uoT  '  (3.3G) 


InterpTetation  of  the  Degrees  of  Freedom 

The  degrees  of  freedom  associated  with  vertices  are  simply  the  function  values  evaluated 
at  these  points.  The  degrees  of  freedom  zissociated  with  the  midpoints  of  element  edges 
and  central  nodes  are  more  complicated.  It  follows  from  definition  (4.36)  that  they  may  be 
interpreted  as  certain  linear  combinations  of  directional  derivatives.  The  form  of  directional 
or  mixed  derivatives  appropriate  for  the  master  element  is  preserved  only  if  the  map  T  is 
linear.  Thus,  it  is  understood  that  the  degrees  of  freedom  discussed  in  the  previous  section 
should  be  interpreted  broadly  enough  to  include  linear  combinations  of  these  Ilermite-type 
degrees  of  freedom. 

Continuity  for  Regular  Meshes 

One  of  the  fundamental  advantages  of  using  the  hierarchical  shape  functions  is  the  ease 
with  which  they  allow  one  to  construct  a  continuous  approximation  with  locally  variable 
order  of  approximation.  A  typical  situation  is  illustrated  in  Fig.  6.  If  elements  A'l  and 
1(2  are  to  support  polynomials  of  degree,  say,  one  and  three,  respectively,  then  there  are  at 
least  two  ways  to  enforce  continuity  across  the  interelement  boundary.  One  way  is  to  add 
two  extra  shape  functions  of  second  and  third  order  corresponding  to  the  middle  node  A  of 
element  K\.  Alternatively,  the  same  two  shape  functions  may  be  deleted  from  element  A'2. 
One  can,  of  course,  establish  a  “golden  rule”  method  by  adding  one  extra  shape  function 
of  second  order  to  Ki,  and  deleting  the  third  order  shape  function  from  A'2.  In  all  these 
cases,  a  common  order  of  approximation  along  the  interelement  boundary  can  be  enforced 
by  simply  adding  or  deleting  the  respective  shape  functions  from  the  neighboring  elements. 
While  any  of  these  choices  can  be  made,  the  results  described  here  employ  the  “maximum 
rule”  in  which  the  higher-order  approximation  dominates  lower  orders.  Thus,  if  an  element 
is  p-refined,  i.e.,  a  higher  order  approximation  of  degree  p  is  introduced,  the  neighbors  of 
lower  order  are  enriched  by  the  addition  of  extra  shape  functions  of  degree  p  necessary  to 
guarantee  continuity  of  the  approximation. 
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Figure  6:  Continuity  by  hierarchical  shape  functions. 
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Constraints  in  the  One- Dimensional  Case 

In  the  case  of  irregular  meshes,  continuity  has  to  be  enforced  by  means  of  the  constrained 
approximation.  To  fix  ideas,  consider  the  generic,  one-dimensional  case  shown  in  Fig.  7. 
The  approximation  on  the  small  elements  [—1,0]  and  [0, 1]  must  match  the  approximation 
on  the  large  element  [—1,1]. 

We  first  choose  the  scaling  factors  Ap  in  (4.25)  in  such  a  way  that  the  corresponding 
shape  functions  for  the  one-dimensional  master  element  have  the  following  form 

Xo  =  ^(1-0 

Xi  =  *  (3-37) 

f  V’-l  P  =  2,4,6,... 

Xp(0  = 

[  e-i  p  =  3,5,7,... 

Assume  next  that  all  degrees  of  freedom  for  the  large  element  are  active.  The  question 
is:  what  degrees  of  freedom  must  the  small  elements  take  on  in  order  that  the  functions 
supported  on  the  two  small  elements  exactly  coincide  with  shape  functions  of  the  large 
element? 

From  the  fact  that  (4.25)  is  a  dual  basis  to  (4.37),  we  get 

V’pCXp)  =  v-p!  =  1  P  =  2, 3, . . .  (3.38) 

and  therefore  Ap  =  p!. 

The  transformation  map  from  [—1,1]  onto  [—1,0]  is  of  the  form 

x  =  -i  +  5e  (3.39) 

with  inverse  ^  =  2x  +  1.  This  yields  the  following  formulas  for  the  shape  functions  ^Xp,p  = 
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Figure  7:  Derivation  of  the  constraint  coefficients. 


0,1,2,...  For  the  (left-hand  side)  element  [—1,0]  (recall  definition  (4.35)). 


^Xo(^)  =  -X 

^Xi(x)  =  I  -t-  1 

^Xp(a:)  =  l-(2x-f-l)P  p  =  2,4,6,... 

^Xp{x)  =  (2x -h  1)'’ -  (2i -f- 1)  p  =  3,5,7,... 

and  the  corresponding  formulas  for  the  degrees  of  freedom  are  (recall  (4.36)). 

<‘ipQ,u>  =  u(-l) 


(3.40) 


<'(pi,u>  =  u(0) 

1  (Pu  (  1' 


,  1  /  1\ 


(3.41) 


Now  let  u(x)  (x  =  ^  for  the  master  element)  be  any  function  defined  by  the  shape 
functions  on  [—1, 1],  i.e.. 


j=0 


(3.42) 


In  order  to  represent  u(x)  for  x  E  [—1,0]  in  terms  of  the  shape  functions  on  [—1,0],  we 
have  to  calculate  the  values  of  the  degrees  of  freedom  (4.41).  We  get 

k 

.  <^ipo,u>  =  ^o,Xo> ^o,X<i  > 

g=l 


=  <  (fo,U> 

<‘(pi,u>  =  ipoiu)  <‘ (pi,Xo> +<Pi{u)  <‘ (pi,Xi  > 
k 

g=2 

1  ^ 

=  «  <  V’l,"  >  +II  '-Rgl  <  > 

^  9=2 


(3.43) 


where 


'/?,1  =  <‘<fl,Xq> 

-I 


1  if  ^  is  even 
0  otherwise 


(3.44) 
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For  p  >  2, 


where 


<‘^Pp,U>  =  ^o{u)  <'  ^p,XO>  +'^^q{u)  <‘  Vp,Xq  > 


<?=1 


=  0  +  EX  <(Pg,U> 

9=1 


Rqp  ^piXq  ^ 

r  0 


=  < 


for  q  <  p 


(-l)P+^  q! 


(3.45) 


for  q  >  p 


2i  p\{q  —  p)\ 

The  same  procedure  applied  to  the  right-hand  side  element  [0,1]  yields  the  following: 


The  transformation  from  [—1,1]  onto  [0,1] 


with  inverse  ^  =  2x  —  1. 


.  =  i  +  ie 

2  2^ 


The  shape  functions  ''XpiP  =  0, 1,2,. . 


’’XoCa:)  =  1  -  X 

’■Xi(x)  =  X 


"Xp{x) 


1  —  (2x  —  1)P  p  even 

(2x  —  1)P  —  (2x  —  1)  p  odd 


The  degrees  of  freedom  ^<Pp: 


<  Vo, «  >  =  ^^(0) 

<  Vl5W>  =  '“(1) 


<  Vp5  “  >= 


1 

2Pp! 


d^u  / 1  \ 
dxP  \2/ 


(3.46) 


(3.47) 

(3.48) 

(3.49) 


(3.50) 
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The  constraints 


1  1  * 

<  ^(po,U  >=  -  <  >  +-  < 

^  ^  q=2 


where 


(3.51) 


and  for  p  >2 


where 


^RqO  =  1 

1  if  9  is  even 

0  otherwise 

(3.52) 

<  >  =  <  <Pl,U  > 

(.3..53) 

k 

<  ‘Ppi'a  >  =  <  Rqp  < 

9=2 

(3.54) 

r  p  _  r 

riqp  — 

’V’p,Xq  > 

0  for 

q  <p 

(3.55) 

,  ^  (p) 

q>p 

(3.56) 

Arrays  ^Rgp  and  q,p  =  0, ...  ,5  are  presented  in  Fig.  8. 

Constraints  for  2-D  Subparametric  Elements 

Since  the  shape  functions  for  the  2-D  master  element  are  defined  as  tensor  products  of 
the  1-D  functions,  the  results  for  the  1-D  case  hold  exactly  in  the  same  form  in  the  2- 
D  situation,  the  only  difference  being  that  the  calculated  constraint  equations  have  to  be 
applied  to  the  proper  degrees  of  freedom  (see  Fig.  9).  It  follows  from  the  rleHnition  of  the 
subparametric  elements  that  the  constraints  coefficients  are  exactly  the  same,  even  when 
the  elements  have  curved  boundaries.  This  follows  from  the  fact  that  the  shape  functions’ 
behavior  in  a  subparametric  element  on  a  part  of  its  boundary  depends  exclusively  upon  the 
deformation  of  the  part  of  the  boundary,  and  therefore,  any  relation  defined  for  the  shape 
functions  in  the  generic  situation  carries  over  immediately  to  the  case  of  two  small  elements 
sharing  an  edge  with  a  large  element,  so  long  as  the  deformation  of  the  edge  is  identical  in 
all  three  elements.  The  situation  is  illustrated  in  Fig.  9. 
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Figure  8:  The  constraint  coefficients  for  a  sixth  order  approximation.  The  unfilled  coefficients 
are  zero. 
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4  Some  Details  Concerning  the  Data  Structure 


In  the  classical  FEM,  elements  as  well  as  nodes  are  usually  numbered  consecutively  in  an 
attempt  to  produce  a  minimal  band  within  the  global  stiffness  matrix.  When  the  program 
identifies  an  element  to  process  its  contribution  to  the  global  matrices,  the  minimal  infor¬ 
mation  needed  is  the  node  numbers  associated  with  the  element.  Adaptive  refinement  and 
unrefinement  algorithms  require  much  more  information  on  the  mesh  structure  than  the 
classical  assembly  process. 

First  of  all,  we  introduce  the  notion  of  a  family.  Whenever  an  element  is  refined  a  new 
family  is  created.  The  original  element  is  called  the  father  of  the  family  and  the  four  new 
elements  are  called  its  sons.  Graphically,  the  geneology  on  families  can  be  presented  in  a 
family  tree  structure  as  illustrated  in  Fig.  10. 

An  examination  of  refinement  and  unrefinement  algorithms  (see  [7]  for  details)  reveals 
that  for  a  given  element  NEL,  one  must  have  access  the  following  information: 

•  element  node  numbers 

•  element  neighbors 

•  the  tree  structure  information,  including:  the 

—  number  of  the  element  family 

—  number  of  the  father 

—  numbers  of  sons 

—  refinement  level  (number  of  generation) 

For  a  given  NODE  we  also  require, 

•  node  coordinates 

•  values  of  the  degrees  of  freedom  associated  with  the  node 

In  general,  some  information  is  stored  explicitly  in  a  data  base  consisting  of  a  number  of 
arrays,  some  other  information  is  recovered  from  the  data  base  by  means  of  simple  algorithms. 
A  careful  balance  should  be  maintained  between  the  amount  of  information  stored  (storage 
requirements)  and  recovered  (time). 

The  following  is  s  short  list  of  arrays  used  in  the  data  base: 

1.  The  tree  structure  is  stored  in  the  condensed,  family-like  fashion  [6],  [7]  in  two  arrays 
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0 
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Figure  10:  A  tree  structure  and  the  natural  order  of  elements:  4,  5,  12,  13,  14,  16,  17,  18, 
19,  7,  8,  9,  10,  20,  21,  22,  23,  3. 
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NSON(NRELEI) 

NTREE(5,MAXNRFAM) 

where  NRELEI  is  the  number  of  elements  in  the  initial  mesh  and  MAXNRFAM  is  the 
anticipated  maximum  number  of  families.  For  an  element  NEL  of  the  initial  mesh, 
NSON(NEL)  contains  its  first  son  number  (if  there  is  any).  For  a  family  NFAM, 
NTREE(1,NFAM)  contains  the  number  of  the  father  of  the  family  while  the  other  four 
entries  NTREE(2:5,NFAM)  are  reserved  for  the  “first-bom”  sons  of  the  sons  of  the 
family  (the  first-born  “grandsons”  of  the  father). 

2.  The  initial  mesh  neighbor  information  is  stored  explicitly  in  array 

NEIG(4, NRELEI) 

containing  up  to  four  neighbors  for  each  element  of  the  initial  mesh  (elements  adjacent 
to  the  boundary  may  have  less  neighbors). 

3.  For  every  active  element,  up  to  nine  nicknames  are  stored  in  array 

N0DES(9,MAXNRELEM) 

where  MAXNRELEM  is  the  anticipated  maximum  number  of  elements. 

For  a  regular  node,  the  nickname  is  defined  as 

NODEnoO  +  NORDER 

where  NODE  is  the  node  number  and  NORDER  the  order  of  approximation  associated 
with  the  node. 

For  an  irregular  node,  the  nickname  is  defined  as 
-  NORDER 

where  NORDER  is  again  the  order  of  approximation  corresponding  to  the  node. 

4.  For  a  particular  component  lEL  of  a  vector-valued  solution,  the  corresponding  degrees 
of  freedom  are  stored  sequentially  in  array 

U(MAXNRDOF,IEL) 

where  MAXNRDOF  is  the  anticipated  maximal  number  of  degrees  of  freedom.  Two 
extra  integer  arrays  are  introduced  to  handle  the  information  stored  in  array  U.  Array 

NADRES(MAXNRNODE) 
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contains  for  every  node,  NODE,  the  address  of  the  first  from  the  degrees  of  freedom 
corresponding  to  NODE  in  array  U.  If  K  =  NADRES(NODE)  is  such  an  address,  the 
address  for  the  next  degree  of  freedom  can  be  found  in 

NU(K) 

and  so  on,  until  NU(K)=0,  which  means  that  the  last  degree  of  freedom  for  a  node  has 
been  found  (see  [7]  for  a  detailed  discussion).  The  parameter  MAXNRNODE  above  is 
the  anticipated  maximal  number  of  nodes. 

5.  The  node  coordinates  are  stored  in  array 
XNODE 

XN0DE(2,MAXNRN0DE) 

The  rest  of  the  necessary  information  is  reconstructed  from  the  data  structure  by  means  of 
simple  algorithms.  These  include: 

—  calculation  of  up  to  eight  neighbors  for  an  element 

-  calculation  of  local  coordinates  of  nine  nodes  for  an  element  determining  its  ge¬ 
ometry  (the  irregular  nodes  coordinates  have  to  be  reconstructed  by  interpolating 
regular  nodes  coordinates) 

-  recovery  of  the  tree-structure  related  information,  e.g.,  level  of  refinement,  the 
sons  numbers,  etc. 

—  an  algorithm  establishing  the  natural  order  of  elements 

During  the  h  and  p  refinements  and  unrefinements,  both  elements  and  nodes  are  created 
and  deleted  in  a  rather  random  way.  This  makes  it  impossible  to  denumerate  them  in  a 
consecutive  way,  according  to  their  numbers  (for  instance,  as  a  result  of  unrefinements  some 
numbers  may  be  simply  missing!).  Thus  a  new  ordering  of  elements  hcis  to  be  introduced 
which  is  based  on  some  scheme  other  than  an  element  numbers  criterion.  In  the  code 
discussed  here,  we  use  “the  natural  order  of  elements”  based  on  the  initial  mesh  elements 
ordering  and  the  tree  structure.  The  concept  is  illustrated  in  Fig.  10.  One  has  to  basically 
follow  the  tree  of  elements  obeying  the  order  of  elements  in  the  initial  mesh  and  the  order 
of  sons  in  a  family. 

The  natural  order  of  elements  may  serve  as  a  basis  for  defining  an  order  for  nodes  and, 
consequently,  for  degrees  of  freedom,  when  necessary. 

For  a  detailed  discussion  of  the  data  structure  as  well  as  a  critical  review  of  different  data 
structures  in  context  of  different  h-refinement  techniques,  we  refer  again  to  [7]. 
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APPENDIX  B 


Theory  and  Implementation  of  High— Order  Adaptive 
/ip-Methods  for  the  Analysis  of  Incompressible  Flows 


1  Introduction 

This  appendix  is  presented  as  a  summary  of  some  approaches  and  some  early  numerical 
results  on  k-p  methods  for  incompressible  flows.  The  appendix  is  based  on  a  forthcoming 
paper  by  J.  T.  Oden,  with  the  same  title  as  that  above,  that  is  to  appear  in  Computational 
Nonlinear  Mechanics  in  Aerospace  Engineering,  a  volume  in  the  “AIAA  Progress  in 
Aeronautics  and  Astronautics.” 

It  is  worth  re-emphasizing  the  strategy  behind  the  uses  of  hp-methodology.  In  impor¬ 
tant  CFD  circles,  the  inadequacy  of  contemporary  computing  devices  for  handling  realistic 
three-dimensional  flow  simulations  is  frequently  presented  as  justification  for  research  on 
and  development  of  larger  and  faster  computations  and  concomitant  larger  budgets  for  com¬ 
puting.  While  few  will  argue  against  the  need  for  bigger  and  better  hardware,  and  we  shall 
certainly  not  do  that  here,  the  premise  for  many  of  these  widely  published  projections  is 
rooted  in  the  jargon  of  traditional,  lower-order,  finite  difference  approaches  for  structured 
meshes.  A  typical  scenario  is  this: 


Problem  Class 

•  3D  Steady  Viscous  Flow 

•  Reynolds  Averaged  Navier-Stokes 

•  Large  Eddy  Simulations 

•  Full  Navier-Stokes 


Number  of  Grid 
Points  Required 

10®  -  10® 

10^ 

10^ 

10«  _  10»5 


Computer  Class 


VI 

30  X  VI 
3  X  10^  X  VI 
3  X  10^  X  VI 


Here  “Class  VI”  refers  to  a  sixth-generation  supercomputer,  such  as  a  CRAY  YMP,  and 
“3  X  10®  X  V/”  suggests  that  machines  with  a  billion  times  the  speed  and  storage  of  today’s 
supercomputers  will  be  needed  to  resolve  flow  features  of  interest  in  high  Reynolds-number 
flows. 

Of  course,  these  figures  are  very  rough  and  one  can  argue  over  the  location  of  decimal 
points  here  or  there,  but  they  do  indeed  highlight  the  enormous  challenge  facing  future  ad¬ 
vances  of  computational  fluid  dynamics.  This  challenge  is  real  and  is  not  at  all  at  issue;  what 
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is  at  issue  is  the  notion  that  flow  features  must  be  resolved  by  merely  adding  and  handling 
more  gridpoints.  This  terminology  suggests  that  the  driving  factor  governing  research  on 
CFD  is  only  the  performance  of  the  computing  device  itself  and  that  the  algorithms  to  be 
implemented  on  these  devices  will  continue  to  be  versions  of  the  fixed  gridpoint  methods 
that  have  been  popular  in  CFD  for  decades. 

An  alternative  philosophy  is  to  also  look  toward  the  development  of  new  algorithms  for 
CFD  that  are  designed  to  streamline  and  optimize  the  computational  process.  These  include 
“smart  algorithms,”  which  change  in  structure  and  performance  during  a  flow  calculation  to 
accommodate  changing  properties  of  the  solution.  In  particular,  they  include  adaptive  finite 
element  methods,  which  are  designed  to  adjust  mesh  parameters  so  as  to  control  numerical 
error.  Importantly,  they  include  so-called  h,  p,  and  r-adaptive  schemes  which  dynamically 
change  the  mesh  size  h,  the  local  spectral  order  p,  and  the  relocation  r  of  gridpoints  so  as 
to  deliver  very  high  accuracies  with  near-minimal  numbers  of  degrees-of-freedom. 

Such  adaptive  methods  need  not  function  on  structured  meshes.  Most  importantly,  they 
can,  and  generally  do,  produce  exponentially-convergent  approximations,  meaning  that  such 
adaptive  methods  can  produce  results  of  a  given  accuracy  using  many  orders-of-magnitude 
fewer  degrees  of  freedom  than  that  required  by  conventional  approaches.  Thus,  flow  features 
are  not  resolved  by  heuristically  adding  gridpoints  where  the  analyst  thinks  they  may  be 
needed,  but  rather  by  automatically  distributing  element  sizes  and  spectral  orders  in  a  way 
designed  to  produce  results  with  a  preset  level  of  accuracy. 

A  remarkable  by-product  of  these  approaches  is  the  calculation  of  estimates  of  the  local 
(elementwise)  error  in  the  approximation.  For  example,  to  control  the  computational  process, 
the  calculated  flowfield  may  be  used  to  derive  element  residuals,  and  these  residuals  can  be 
used  to  calculate  estimates  of  error  in  appropriate  norms.  The  quality  and  sophistication 
of  existing  estimation  methods  vary,  but  they  have  proven  to  be  invaluable  in  giving  the 
analyst  an  independent  measure  of  the  reliability  of  his  calculation.  Also,  they  provide  a 
rational  alternative  to  judging  the  convergence  of  a  flow  calculation  by  simply  monitoring 
the  change  in  results  due  to  successive  refinements. 

While  recognition  of  the  significant  potential  of  adaptive  ftp-methods  is  increasing  in  the 
CFD  community,  many  research  issues  remain  to  be  resolved.  The  present  paper  outlines  a 
class  of  adaptive  ftp-methods  for  incompressible  viscous  flow  problems  and  discusses  results 
on  the  development  of  ftp-strategies,  error  estimators,  higher-order  flow  solvers,  and  results  of 
representative  numerical  experiments.  This  work  continues  our  study  of  adaptive  techniques 
initiated  in  the  early  1980s  for  both  compressible  and  incompressible  flow  simulations. 

Work  on  ft/>-adaptive  methods  in  computational  fluid  dynamics  began  in  1983  and  1984 
(see  [1,2])  and  originally  focused  on  two-dimensional  steady,  viscous  incompressible  flow. 
Improvements  in  techniques  for  error  estimation  and  in  ftp-data  structures  led  to  extensions 
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to  compressible  flow  problems  (e.g.,  [3-8]);  for  surveys  on  general  adaptive  finite  element 
methods,  see  [9-11].  Smart  algorithms  and  adaptive  methods  in  CFD  are  discussed  in  [12- 
16].  An  overview  of  h-  and  h-p  finite  element  methods  for  elliptic  problems  is  given  in 
[17]. 


/ip-finite  element  methods  share  features  of  standard  fixed-order  finite  element  methods 
and  of  spectral  methods,  both  of  which  have  a  long  history.  For  a  history  of  finite  element 
methods,  see  [18].  Global  spectral  methods,  of  course,  data  back  to  the  nineteenth  century 
and  were  used  in  a  weak  or  variational  setting  by  Ritz  [19]  in  1908,  Galerkin  [20]  in  1915,  and 
Kantorovich  and  Krylov  [21]  in  1936.  The  idea  of  using  spectral  approximations  for  finite 
element  methods,  with  or  without  collocation,  appeared  very  early  in  the  finite  element 
literature  (see,  e.g.,  [21]).  The  idea  of  using  exact  or  high-order  representations  of  functions 
over  a  regular  subdomain  and  to  then  patch  these  subdomains  together  to  produce  a  global 
model  was  the  forerunner  to  “more  practical”  finite  element  methods  based  on  low-order 
approximations  and  finer  meshes.  The  problem  was  with  interface  conditions — how  could 
these  high-order  spectral  representations  be  effectively  matched  together? 


The  interface  issue  was  rather  trivially  solved  when  only  low-order  conforming  polyno¬ 
mial  approximations  were  employed,  and  this  fact  led  to  the  development  of  traditional  and 
highly  successful  finite  element  schemes.  The  resurrection  of  high-order  schemes  within  a 
finite  element  context  came  with  the  pioneering  work  of  Szabo  and  his  collaborators  (e.g. 
[23-28])  on  the  p- version  of  the  finite  element  method.  There  the  notion  of  hierarchical  fam¬ 
ilies  of  polynomials  on  triangular  meshes  was  developed,  which  simultaneously  resolved  the 
interface  problem  and  led  to  techniques  which  produced  remarkable  accuracy  for  classes  of 
elliptic  problems  in  linear  elasticity.  The  mathematical  properties  of  p- version  finite  elements 
has  been  studied  by  Babuska,  Szabo,  and  others  [17,24,29-35]  and  is  now  a  widely-accepted 
approach  for  solving  problems  m  linear  elasticity  and  in  other  areas.  The  spectral  element 
methods  developed  by  Patera  and  Maday  and  their  colleagues  (e.g.,  [36-40])  are,  in  fact, 
closely  related  to  p-versions  of  the  finite  element  technique.  In  this  technique,  Chebyschev 
polynomials  cind  either  collocation  or  Galerkin/Gauss-Lobbatto  quadratures  are  used  to 
produce  high-order  local  spectral  approximations  on  a  fixed  mesh.  Interface  conditions  are 
enforced  by  point  matching  at  boundary  collocation  points  or  by  using  a  “weak”  continu¬ 
ity  condition  in  which  interelement  continuity  is  enforced  by  an  L^-projection  of  traces  of 
functions  on  neighboring  elements  onto  their  common  boundary. 


A  general  approach  to  /ip-adaptive  methods  was  developed  in  a  series  of  papers  by  Oden, 
Demkowicz,  Rachowicz,  and  others,  see  [41-43].  There  a  general  data  structure,  constructed 
on  quadrilateral  or  hexahedral  meshes,  was  developed  which  provided  for  irregular  h  refine¬ 
ments  and  p  enrichments.  An  adaptive  strategy  is  also  given  in  these  works  together  with  a 
mesh  optimization  scheme  designed  to  choose  an  acceptable  sequence  oi  h  or  p  refinements 
to  yield  an  optimal  nonuniform  hp  mesh.  An  extension  of  several  aspects  of  these  works  is 
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given  in  the  present  paper. 

Following  this  Introduction,  mathematical  preliminaries,  notations,  and  weak  forms  of  the 
Navier-Stokes  equations  are  given  in  Section  2.  General  properties  of  /ip-finite  element  spaces 
are  established  in  Section  3,  and  /ip-finite  element  models  of  the  Navier-Stokes  equations 
are  developed  in  Section  4.  There  we  also  discuss  penalty  approximations  of  the  pressure 
and  incompressibility  constraint,  which  are  used  in  developing  higher-order  flow  solvers. 
Several  high-order  semi-implicit  solvers  are  discussed  in  Section  5  and  error  estimation  and 
adaptivity  techniques  are  discussed  in  Section  6.  Section  7  contains  the  results  of  several 
numerical  experiments. 


2  PRELIMINARIES 


In  this  section,  we  record  preliminary  notations  and  problem  definition  for  future  refer¬ 
ence. 

2.1  The  Governing  Equations 

We  begin  with  the  description  of  the  flow  of  a  viscous  incompressible  fluid  embodied  in  the 
following  problem; 
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Given  an  open  bounded  domain  Q  C  i8",  n  =  2  or  3,  data 
(f  o,nd  time  interval  [0,7’],  find  a  velocity  field  u  = 

u{x,t)  and  a  pressure  P  =  P{x,t),  x  =  (xi,x2,  •  •  •  ,x„)  6  H, 
t  G  [0,  T],  such  that 


Ut  +  u-  Vu  -  2i/ V  •  D{u)  +  VP  =  f 

V  •  tt  =  0 


in  n  X  (0,  T] 


M(a;,0)  =  ■uo(®)  in 
u{s,t)  =  u{s,t), 

se^^,  <€[o,r] 

o-(u,P)-n(s,t)  =  g(s,t), 

5 era,  <€  [o,r] 


Here  the  following  notations  and  conventions  are  used: 


du 

Ut 

dt 

D{u) 

= 

—  (V'u  +  =  the  strain  rate  tensor 

V 

= 

hIp  =.  the  kinematic  viscousity,  y,  being  the  constant  viscosity 

• 

and  p  the  density  of  the  fluid 

p 

= 

p! p  =  the  kinematic  pressure;  p  being  the  hydrostatic  pressure 

f 

— 

f{x,t)  =  fjp  =  the  body  force  per  unit  volume,  /  being  the 
body  force  density 

Uo 

= 

Uo{x)  =  the  initial  velocity  distribution 

u 

= 

«(s,<)  =  prescribed  velocity  on  a  portion  Fj  of  the  boundary 
of  fl: 

X  6  do,  =r'  X  =  ®(a),  Fi  C 

tr{u,P) 

= 

the  (density-weighted)  Cauchy  stress  tensor  =  p{2i>D[u)  — P\), 
1  being  the  unit  tensor 

n 

= 

a  unit  exterior  normal  to  dfi 

9 

= 

g(s,t)  =  prescribed  traction  on  F2  C  dfl, 
dQ  =  riU  F2,  Fi  n  Fi  =  0 

Remarks: 

2.1)  The  boundary  and  initial  data  must  satisfy  certain  compatibility  conditions  if  (2.1) 
is  to  represent  a  well-posed  problem.  In  particular,  we  must  have 


V  •  tio  =  0  and  «o  =  wC-SjO)  on  Tj 


Uo  ■  n  ds  =  0 


We  demand  that  meas  Fj  >  0  or  F,  =  0,  i  =  1,2,  then,  if  F2  =  0,  w  must  satisfy. 


I 


u  ■  n  ds  —  0 


ri=an 


te  [0,T] 


an 


2.2)  Boundary  conditions  other  than  those  given  in  (2.1)  can  be  incorporated  into  the 
alysis.  See,  e.g.,  [44].  | 


2.2  Weak  Formulations 

To  recast  (2.1)  into  a  functional  setting,  we  introduce  the  following  notations: 
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Spaces 


V  =  |i;  =  t;(a)  G  (/f*(n))”  I70W  =  0  on  Fij 
H  =  {r  G  V\div  v  =  0  in  fl} 

M  =  {qe  L\n)\  f  q  dx  =  oifr2  =  0] 


(2.2) 


When  the  domain  of  functions  needs  emphasis,  we  will  also  use  the  notations  V  =  V{n), 
H  =  H{U),  and  M  =  M(fl).  We  also  introduce  the  norms, 

=  J^CVv  :  "Vv +  v  •v)dx, 
ll9llo,n  =  L  ll^llr.n  =  >  ^  >  ^tc. 

t=i 


Forms 


a:  VxV-^M 

a{u,v)  =  f  2iyD{u)  :  D{v  )  dx 


b:  VxM-*R 

b{v,q)  =  q  div  v  dx 
Jo 

c:  VxVxV 

c(u,  V  ;  w)  =  J  w  dx 

L  :  V^E 

L(^)  =  (/>'*')  +  /  g-v  ds  J 


(2.3) 


Here  /f^(fi)  is  the  usual  Sobolev  space  of  functions  with  L^-distributional  derivatives  defined 
on  fi,  7o  is  the  trace  map  of  (n))”-functions  onto  dfl  (jqV  =  v\d^,  v  G  (C°  (f^))  , 
and  we  shall  simply  write  701;  =  t;  unless  confusion  is  likely),  and  (/,  v)  =  Sq  f  ■  v  dx, 
dx  =  dx\dx2  ■  •  •  dxn-  Here  also  the  notations  A  :  B  =  Yl7,j=i  and  a  ■  A  ■  b  =  Y17.j=i 

CiAijbj  are  used  for  tensors  A,B  and  vectors  a,  6.  A  weak  formulation  of  (2.1)  can  now  be 
presented  as  follows: 
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Given  fe  L^{0,T-,L\n)),g  €  L\T2)),  uo  G  H, 

and  it  G  (77^(0))",  u|ri  =  «,  find  u  =  u  +  L'^{0,T;V) 
and  P  G  T ;  M),  tt(0)  =  Mq?  such  that 


{ut,v)  +  c{u,u;v)  +  a{u,v)  —  b{v,P)  =  L{v) 

b{u,q)  =  0 

V(t?,9)eV  X  M 


(2.4) 


Conditions  on  the  regularity  of  the  data  {f,g,Uo,  u)  can,  of  course,  be  weakened;  see  Temam 
[45]. 

The  Stokes  problem  arises  as  a  special  case  of  (2.4)  when  Ut  =  0  and  the  convection  term 
c('  ,  •  ;  •)  is  dropped: 


2.3  Weak  Forms  on  Partitions  of  Q, 

It  is  convenient  to  also  define  weak  formulations  of  the  Navier-Stokes  equations  on  partitions 
of  n.  Let  Qh  denote  a  family  of  partitions  7^  of  into  subdomains  K: 

U=\J{KeT,},  Kf]J  =  rKj  v/r,jGT, 

h>0 

hfi  =  dia  {!<),  h  =  sup  A/c,  %  G  Qh 

Ken 
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We  denote  by  V{K),  H{K),  and  M{K)  the  spaces  of  restrictions  of  functions  in  V,H,  and 
M,  respectively,  to  subdomain  K.  Let  (•,•),  ;  Oi  and  Lk{-)  denote  forms 

analogous  to  (2.3)  defined  on  these  restricted  spaces;  i.e., 


a/f(u,  v)  =  /  2i/D{u)  :  D{v)  dx 
J  K 

bj^{v,q)  =  j  q  div  v  dx 
•J  IC 

ck(u,  v;  w)  =  I  (u  ■  Vi>)  •  w  dx 
J  K 

LK{v)=={f,v)ic+  I  gvds 
•/rjnsic 


with  (/, v)k  =  fKf-v  dx.  Then  a(u, v)  =  Eifer^  o/f («,  v),  b  (v, q)  =  ?)’  •  •  •’ 

etc.,  and  the  weak  form  of  the  Navier-Stokes  equations  for  a  typical  subdomain  K  assumes 
the  form. 


We  construct  approximations  of  (2.4)  by  using  /ip-approximations  of  the  spaces  V  and  M  in 
Section  4. 
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3  hp  FINITE  ELEMENT  APPROXIMATIONS 

We  now  take  up  the  fundamental  issue  of  constructing  /ip-finite  element  approximations 
of  the  spaces  V,  H,  and  M  and  V{K),  H{K),  and  M{K). 

3.1  Shape  Functions 

To  simplify  the  discussion,  we  begin  with  the  two-dimensional  ca^e  illustrated  in  Fig.  1. 
There  we  see  the  domain  partitioned  into  a  mesh  of  elements  K  each  of  which  is  the  image 
of  a  master  element  under  a  smooth  invertible  map  Fk'- 

X  =  {xi,X2)  =  Fk{^)  ,  ^  =  (6,6) 


Let 

^  =  [-l,l]x[-l,l] 
We  next  record  properties  of  special  shape  functions  on 
Lemma  3.1.  Let  ^  €  [—1, 1]  and  denote 


1  ~  € 

=  m  f  Pj^i{s)ds  J  >  2 


(F/(0  - 


\/2(2j-l) 

where  Pj{^)  is  the  Legendre  polynomial  of  degree  j.  Then 

G  ■Pj([— 1, 1])  =  the  space  of  polynomials  of  degree  j  ora  [—1, 1]  'j 

<p,(±l)  =  0 

r  [{2j  +  1)-1  +  (2j  -  3)-‘]  /2(2j  -  1)  if  j  =  k 


{T3.Tk)=  < 


-  (2;  +  1)-^  /2  [(2j  -  l)(2j  +  2)]^  if  j  =  k-2 


-  {2j  -  3)-'  /2  [i2j  -  l){2j  -  if  j  =  k  +  2 


if  otherwise 


[{2j  -  l)/2]/{2j  -  1)  if  j  =  k 


if  otherwise 


(3.1) 


(3.2) 
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where  (9,^’)  = 

Proof:  These  results  follow  from  well-known  properties  of  Legendre  polynomials.  See  [46] 
and,  particularly,  [47].  B 

With  these  notations  in  hand,  we  define  the  master  element  shape  functions  in  three 
groups: 

Node  Functions.  These  are  the  standard  bilinear  shape  functions 


=  A=1,2,3,4  (3.3) 

U  =  =  6) 

Side  Functions.  On  each  side  Fa  of  the  boundary  dK  of  /v^,  we  define 

=  5(1 +  >!)«({) 

xigU.v)  =  5(1 -OwW 

xl'~((,v)  =  5(i--i)w(0  (S'*) 

x'ggdiv)  =  j(i+Ov’»(’;) 

k  =  2,3,*-*,p 

Bubble  Functions.  Internally,  we  define  functions  which  vanish  on  dK  by 


xfkii^'n)  =  ,  2<j,k<p  (3.5) 


The  space  \  l{j  of  master  element  shape  functions  now  defined  by 

(k)  =  span  X%  •  1  <  F,  A  <  4,  2  <  j,  k  <  p}  (3.6) 


60 


In  three  dimensions,  similar  choices  of  shape  functions  are  used  except  that  there  it  is 
convenient  to  decompose  the  set  into  edge  and  face  functions  tis  well  (see  Fig.  2).  For 
instance,  for  K  =  [-1, 1]^,  ^  =  6>  »7  =  6,  C  =  6.  we  have 

3D  Node  Functions 


7, 0  =  ^(1  ±  0(1  ±  0(1  ±  0  1  A  =  1, 2,  •  ■  ■ ,  8 


(3.7) 


Edge  Functions 


Face  Functions 

x]k(tv,0 

x%{LvX) 


j(l  +  {)(1  +  CW-!) 
j(i  +0('  -Ovtii) 

fc  =  2, 3,  •  •  • ,  p 

^(i  +  0xf*(7,C) 

^(1  -Oxf*(»/,C) 

^(1  -0x®(00 

‘2-<  j,k  <p 


Bubble  Functions 

xfkeU,V,0  =  <Pj{0'Pk{r})M0 

2  <  j,  <  p 


(3.S) 


(3.9) 


(3.10) 
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Then  we  denote  the  space  of  master  element  shape  functions  by 


(A')  =  span  Xe  ^  ^  <  8, 1  <  A  <  12, 1  <  F  <  6, 2  <  >,  -t,  ^  <  p} 

Denote  dim  Q^{K)  =  N  =  (p  +  1)^.  Then  each  2  €  Q^(K)  representable  in  the  form 

“(0  =  5Za'0.(O 

t=i 

with  0,  the  (relabeled)  shape  functions  defining  the  basis  of  Q^{K). 

3.2  Subparametric  Maps 

Returning  to  the  case  (3.3)  for  convenience,  we  choose  for  the  element  maps  Fk  the  sub- 
parametric  functions, 

*  =  ^K'(^)  =  FE{(,Tf) 

4  4  s  (3.12) 

=  aK+ j^XAX^(6»?)  +  EE®AXa^(^,»?) 

Aal  r=l k=2 

where  is  a  constant  translation  vector,  Xa  =  (xf ,  xf)  are  the  coordinates  of  node  (vertex) 
A  of  element  K,  ajid  Xa  =  (2:^1  ^re  the  coordinates  of  points  interior  to  the  edges 
of  dK .  Often  we  take  only  5  =  2  so  that  element  geometry  is  modeled  using  subparametric 
quadratic  approximation  of  curvilinear  boundaries. 


3.3  h-p  Meshes  and  Data  Structure 

The  element  shape  functions  described  above  are  hierarchical,  meaning  that  the  parameter  p 
can  be  increased  or  decreased  with  the  result  that  the  Jissociated  element  matrices  for  p  =  Pi 
are  always  submatrices  for  p  =  p2,  P2  >  Pi-  We  also  treat  the  mesh  size  h  {hx  =  dia{K)) 
as  a  parameter  by  employing  standard  A-adaptive  techniques  for  bisection  of  elements  dis¬ 
cussed,  for  example,  in  [4].  Conformity  of  shape  functions  across  interelement  boundaries 
is  accomplished  using  constraint  approximations  of  the  type  introduced  in  [41].  The  basic 
ideas  are  summarized  as  follows: 
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1.  Let  7^  G  be  a  partition  of  into  M  =  M{h)  elements  /\, =  U{/\  6  7^}.  Let 
S^^{K)  denote  the  space  of  local  finite  element  functions  defined  over  element  /v , 

l  t=i 


X  e  K  ,  ifi  =  ^iO  ,  ipi  e  (/i) , 

em  I 


where  Q^{K)  is  defined  in  (3.6). 

2.  Let  S(/f)  =  denote  the  space  of  degrees-of-freedom  corresponding  to 

S’^^{K);  i.e.,  T,{K)  can  be  identified  as  the  dual  space  of  S^’’{K)  so  that  II(A')  contains 
hnear  functionals  q)^.  In  particular,  the  functionals  q)^  can  be  chosen  to  be  a  dual  basis 
to  the  basis  functions  V’t  €  S^^{K): 

(qk^'Pj)  =^j  y  1  <  ^ 

where  (•,  •)  denotes  duality  pairing  on  E(A’)  x  S^^{K).  Then,  if 


t;  €  S^^{K) 
v{x)  =  ^v'xl>i{x) 

t=i 

we  must  have 

=  Qk(v)  =  (9a:,  v) 

The  functionals  q}^  (the  degrees-of-freedom)  consist  of  nodal  values, 

9k(«)  =  «(*k)  a  =  1,2, 3, 4 

and  directional  derivatives  at  the  sides  or  mixed  derivatives  at  the  centroid  of  A; 

e.g., 

q:  u — >  ,  A:  =  0,1,  ••• 

where  Z?^u  is  the  ^-th  order  differential  of  u. 

3.  The  space  A''(f2)  of  unconstrained  degrees  of  freedom  consists  of  the  set  of  functions 
Uh'.  Vt  R  such  that 

i)  Uh\K&S^^{K)  VKeTn 


64 


and 


ii)  for  every  pair  of  elements  K,J^Th  sharing  a  node  G  K  HJ,  and  every 
pair  of  degrees  of  freedom, 

qK  €  qj  €  S(J)  with 


A:  =  0, 1,  •  •  • ,  we  have 


quiu)  =  qj{u)  D^u  (x^)  •  (^1 ,  •  ■  • . 


4.  The  global  degrees  of  freedom  E(f2)  are  linear  functionals  on  X(n): 

S:  Xin)-^E,  Hiu)  =  qK{u\K) 

and  the  {discontinuous)  global  base  functions  Xk  are  defined  so  that 


{i:\xk)  =  Si,  i<j,k<N 


N  =  dim 


5.  Now  a  conforming  set  of  finite  element  basis  functions  x  a^re  obtained  by  imposing 
continuity  constraints.  For  this  purpose,  we  partition  the  degrees  of  freedom  into  sets 
N'^  and  of  active  and  constrained  (indices  of)  degrees  of  freedom,  respectively. 
Then,  for  each  constrained  degree  of  freedom  S’,  i  6  N^,  there  is  a  set  I(i)  C  such 
that  (see  [41]) 

/€/(!) 

where  R)  defines  the  Bookean  map  of  active  into  constrained  sets.  For  K,  the  entries  i?) 
are  rational  numbers  the  exact  nature  of  which  depends  upon  the  fi-adaptive  strategy. 

The  final  global  finite  element  space  is  then 


=  ( 


X'’P(fl)  =  {u;,GX'‘(O)|S’(v0=  E  ViG  AT- 


X.(®) 


=  span{x,} 

=  X.(®)  +  E  ^iXj  , 

i€5{i) 


(3.13) 


*^(0  =  {;  e  N^\i  €  /(i)} 
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3.4  The  Zip-Finite  Element  Spaces 

The  construction  outlined  in  the  previous  subsection  leads  to  the  definition  of  families  of 
finite-dimensional  spaces  of  functions  of  the  type  (for  the  two-dimensional  case) 

x^p(n)  c  H\n)  ,n  =  \J  {k  e 

=  (v  S  C»(n)  1 €  S‘'( A-)  © 

I  fc=l 

A/‘-(n)  =  {,  €  i^fi)  I  ,1k  e  5*-(/o®  Bjca/C) 

I  it=l 

We  use  the  notations: 

Pk  =  the  order  of  the  highest-degree  complete  polynomial  rep¬ 

resented  in  the  shape  functions  defined  on  element  K 

Ef'{dK)  =  the  spaces  of  functions  of  the  type  (3.4)  defined  on  the 
edges  of  K  which  are  introduced  to  fulfill  continuity 
requirements  across  element  boundaries 


Thus,  for  a  typical  constrained  element  K  in  the  mesh,  such  as  is  shown  in  Fig.  3,  continuity 
of  the  global  basis  functions  is  maintained  by  enriching  the  edge  functions  of  K  to  match 
the  highest  order  polynomial  experienced  on  the  common  interelement  boundary.  When 
adjoining  element  edge  functions  are  polynomials  of  different  degree,  the  edge  functions  on 
common  boundaries  are  enriched  by  the  addition  of  edge  functions  of  a  degree  equal  to 
the  largest  polynomial  degree  of  any  element  sharing  that  edge.  Thus,  for  the  hp-mesh 
shown  in  Fig.  4,  shape  functions  of  a  given  degree  may  overlap  those  of  adjacent  elements 
which  support  lower  degree  shape  functions.  This  is  illustrated  in  the  color-coded  figure  by 
colored  patches  extending  over  common  element  edges.  In  three-dimensional  /ip-meshes,  the 
implementation  of  this  strategy  requires  also  the  addition  of  element  edge  functions  as  well. 
Once  continuity  is  established,  the  highest  order  of  complete  polynomials  contained  in  the 
element  is  denoted  pj^. 

Typical  interpolation  properties  of  ^p-spaces  are  indicated  in  the  following: 

Proposition  3.1.  Consider  a  mesh  Th  of  elements  K  which  are  affine  equivalent  to  a 
master  element  K ,  with 


K  =  [-hlf  ,  n  =  \J{K  ,  Fk-.K^K 


K  3  X  =  Tkx  +  aK  , 


(3.15) 
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Figure  3:  Illustration  of  the  addition  of  edge  polynomials  at  constrained  edges  so  as  to 
enforce  continuity  of  the  basis  functions  across  interelement  boundaries. 
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Figure  4:  Color-coded  nonuniform  /ip-mesh. 
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Tk  being  an  invertible  (constant)  matrix  and  ax  a  translation  vector.  Let 

px  —  sup{dia  5  I  5  =  sphere  contained  in  K) 
hx  =  dia  [K) 

and  suppose 

— ^  <  (T  =  const.  V  K  E  Th 
PK 

Then,  given  any  function  u  €  there  exists  a  sequence  of  interpolants  E  Vp,^{K), 

the  space  of  polynomials  of  degree  <  px  defined  on  K,  px  =  -  ,  and  a  constant  C, 

independent  of  u,px,  or  hx,  such  that  for  any  s,  0  <  s  <  r, 


u-wH  <C%7  ||u|U. 

\\3,k  - 


(3.16) 


where 


p  —  min(pK  +  l,r) 


(3.17) 

Proof:  See,  for  example,  [17].  I 

Returning  to  (2.1),  we  assume  that  0,  is  such  that  families  of  spaces  V^'’(n),  iU^*(H)  can 
be  constructed  such  that 


1. 


c  K(n) 


M'‘"(n) 


c  M{n) 


(3.18) 


2.  V  v  e  V{n)  fl  (^’■(fi))" ,  3  e 

such  that 


V  X  — 


K 


(3.19) 


3.  V  E  Af**(n),  (3.16)  holds  with  px  replaced  by  sx 
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4  FINITE  ELEMENT  FORMULATIONS 


4.1  /ip-Element  Strategies 

A  semi-discrete  /ip-finite  element  approximation  of  the  Navier-Stokes  equations  consists  of 
seeking  velocity  and  pressure  approximations 

for  t  G  [0,  T]  such  that 


(4.1) 


where  and  are  defined  in  (3.18). 

To  proceed  further,  the  following  issues  must  be  resolved: 

i)  Temporal  Approximation'^  the  discretization  of  [0,r]  and  the  method  of  approximating 

the  time- variation  of  and 

ii)  Convection  Terms]  methods  for  handling  the  nonlinear  convection  effects  embodied  in 
the  form  c(-  ,  •  ;  •). 

iii)  Pressure  Approximations]  methods  for  constructing  stable  approximations  on  hp- 
meshes. 

iv)  Diffusion  Terms]  methods  for  handling  the  diffusion  terms,  a(-  ,  •). 

In  the  present  exposition,  these  issues  are  dealt  with  as  follows: 

i)  High-Order,  Semi-Implicit  Splitting  Schemes.  High-order  temporal  approximations 
are  introduced  to  balance  the  high-order  spatial  approximations  possible  in  adaptive  hp- 
methods.  Splitting  methods  are  used  so  that  convection  steps  are  handled  using  explicit  (or 
“weakly”  implicit)  methods,  whilst  diffusion  steps  are  handled  by  a  linear  (parallelizable) 
implicit  step.  These  are  discussed  in  Section  5. 

ii)  Exterior  Penalty  with  Pressure  Correction.  A  non-singular,  perturbed  Lagrangian 
method  is  used  to  approximate  element  pressures  which  is  essentially  a  regularized,  dis¬ 
continuous,  exterior-penalty  type  approximation  of  the  pressure  P.  In  this  way,  the  local 


-fa  (it'‘P(t),x) 

-  =i(x) 

6(tt'‘P(<),(^)  =  0 

Vx6  V''‘P(fI),v?  G  M'**(H) 
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spectral  orders  s  of  the  element  pressures  are  determined  by  the  quadrature  rule  chosen 
to  integrate  the  bilinear  form  b{-  ,  •).  The  regularization  (penalty)  parameter  £  =  de¬ 
pends  upon  /iff  and  pff  and  can  vary  element  to  element.  The  incompressibility  condition, 
div  It  =  0,  is  thus,  relaxed,  with  =  —ej^div  at  quadrature  points  G  K. 

Whenever  Hc/it;  exceeds  a  preset  tolerance,  the  pressure  is  globally  corrected.  These 

details  are  discussed  more  fully  in  the  next  subsection. 

4.2  Pressure  Approximations 

We  begin  with  a  brief  review  of  the  concepts  behind  non-singular  perturbations  in  the 
continuity  equation  as  a  device  for  handling  pressure  approximations. 

It  is  well  known  that  a  notorious  difficulty  in  the  numerical  solution  of  the  Navier-Stokes 
equations  is  the  treatment  of  the  pressure  approximation  in  a  way  that  results  in  a  stable 
scheme.  At  the  heart  of  this  issue  is  the  LBB-condition,  originally  established  for  the  Stokes 
problem  by  Ladyszhenskaya  [48],  which  establishes  that  if  the  Stokes  problem  is  to  be  well 
posed,  the  bilinear  form  b:V  x  M  R  ol  (2.3)  must  be  such  that  a  constant  /3  >  0  exists 
such  that 

sup  V,SM  (4.2) 

V^V\Q 

A  fairly  large  literature  exists  on  ways  to  ensure  the  fulfillment  of  this  condition  for  certain 
discrete  problems  (for  summaries,  see  Carey  and  Oden  [49],  Kikuchi  and  Oden  [50],  and 
Brezzi  and  Fortin  [51]). 

Focusing  on  the  Stokes  problem  (2.5)  for  simplicity,  it  is  observed  that  the  solution  (■u,p) 
is  a  saddle  point  of  the  Lagrangian 

L  :  V  xM-*£ 

L{v,q)  =  J{v)-b{v,q)  >  (4.3) 

J{v)  =  ^a{v,v)- L{v)  \ 

A  unique  saddle  point  of  L  exists  whenever  the  energy  functional  J{v)  is  coercive,  convex, 
and  differentiable  (with  respect  to  v)  for  fixed  q  (which  is  always  the  case  for  the  class  of 
problems  considered  here),  when  L{vo,q)  is  concave  and  differentiable  with  respect  to  q  for 
fixed  V  =  Vo  (which  also  holds),  and  when  the  following  coercivity  condition  holds  for  some 
V  =  Vq  (see  Oden  [52,  p.  81]): 

lim  L{vo,q) — ^ -oo  (4.4) 


71 


To  enforce  this  coercivity  condition,  and  ensure  a  well-posed  problem  for  evaluating  the 
pressure,  one  can  perturb  the  Lagrangian  by  the  addition  of  terms  that  guarantee  that  (4.4) 
holds,  e.g., 

L,{v,  q)  =  L(v,  q)  +  D^{v,  q)  (4.5) 

where  Dc{v,q)  is  designed  so  that  Lj(-  ,  •)  satisfies  (4.4).  Then  there  exists  a  sequence  of 
saddle  points  {ue,Ps)  such  that 


Le  (Ws,  q)  <  (Wo  Pc)  <  ie  (»,  Pc)  V  (r,  9)  e  K  X  M 


and  one  structures  £)j(-  ,  •)  so  that  {ue,Pg)  —*■  (u,P)  as  £  — +  0,  (u,P)  being  a  solution  of 
(2.5). 

For  the  Stokes  problem,  one  of  the  most  common  forms  of  the  perturbation  term  in  (4.5) 
is 

De{v,q)  =  -\\\q\\l  (4.7) 

Then  (4.4)  holds  for  any  £  >  0  and  the  perturbed  form  of  the  Stokes  problem  is 

a{Ug,v)  -  h{v,Pg)  =  L{v)  VreV  1 

(4.8) 

{ePcq)  +  h{uc,q)  =0  VgeMj 

which  is  formally  equivalent  to  an  exterior  penalty  approximation  of  Pj,  since  (4.8)2  yields 


Pg  =  —e~^div  Ug 


With  (4.8),  the  momentum  equation  becomes 


a  {Ug,  v)  e  ^  {div  Ug,  div  v)  =  L{v)  V  r  G  K 


(4.10) 


In  discrete  approximations  of  (4.9),  we  use  numerical  quadrature  rules  for  evaluating  the 
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integral  of  the  penalty  terms.  For  the  master  element  K, 


fgdi  w  /(/,  ^)  =  (^/)  9  ((<) 


f=i 


(4.11) 


where  W/  are  the  quadrature  weights  and  are  the  quadrature  points,  and  for  affine- 
equivalent  elements  (recall  (3.15)),  we  have 


IU,9)=  H  IkU.s)-,  lKU,9)=iciTKHf,9)  (4.12) 

The  space  of  discontinuous  element  pressures  A/^*(fi)  is  defined  by 

A/“({1)  =  {?‘'£M|«*>|^oFk  =  ?6<3‘(a’),  /(?,?)  =  IIC.(k)}  ('*•13) 


Thus,  the  pressures  are  discontinuous  polynomials  of  degree  s  in  each  variable,  with  s  de¬ 
termined  by  the  quadrature  order  of  /(•  ,  •).  For  an  Ap-approximation  of  velocities,  the 
discrete-penalty  formulation  becomes. 


a  +£  (div  u^^,div  =  L{v)  V  r  G  V  (4-14) 


and  the  pressure  approximation  is  given  by 


G  g^n);  n'*"  o  Fk  iCf)  =  -e-Uiv  (^,)  (4.15) 
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The  discrete  LBB  condition  for  this  type  of  approximation  is  of  the  form, 


I  fdiu 

II  hell - - 

t7'>pgv'>p\o  lp"^lli,n 


o.n  ^qeM^^n) 


(4.16) 


A  study  of  LBB  conditions  of  this  type  for  high-order  /ip-finite  element  methods  has 
been  conducted  by  Wu  and  Oden  [53]  for  some  special  two-dimensional  cases.  The  following 
result  suggests  how  a  stable  approximation  of  pressure  can  be  obtained. 

Let  Q  denote  a  rectangular  domain  in  that  is  partitioned  into  a  mesh  of  quadrilateral 
elements  K  which  are  affine  equivalent  to  a  master  element  K  =  [—1,1]^.  Let  the  velocity 
approximation  of  LI  belong  to  the  space, 


=  [o'"  €  6 


(4.17) 


where  Q^{K)  is  the  space  of  hierarchical  polynomial  shape  functions  of  degree  p  in  x\  and 
X2,  and  let  the  pressure  approximation  q^*  belong  to 


G  L^Ll)  G 


(4.18) 


Then  the  discrete  LBB  condition  (4-i5)  is  satisfied  with  0hp  independent  ofh  and p  whenever 


p>2,  s  =  p  —  r,  r  >2 


Numerical  experiments  suggest  that  this  result  also  holds  for  three-dimensional  cases. 


4.3  Pressure  Corrections 

A  great  deal  of  attention  has  been  given  to  pressure-based  schemes  for  treating  the  incom¬ 
pressible  Navier-Stokes  equations  in  recent  years.  These  families  of  methods  originated  in 
the  work  of  Chorin  [54],  were  the  basis  of  the  Simple  and  Simpler  finite  difference  algorithms 
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of  Patankar  and  Spalding  [55]  (see  also  [56]),  and  were  studied  in  great  detail  within  the 
context  of  finite  elements  by  Gresho  (e.g.,  [57,58]);  see  also  Ramaswamy,  et  al.  [59]. 

The  ideas  are  based  on  extensions  of  the  Helmholtz  theorem: 

For  any  vector  field  u:Cl  — +  there  exist  scalar  and  vector  potentials  and  rj) 

such  that 

u  =  +  curl 

Since  curl  =  0  and  div  curl  ij)  =  0,  u  can  be  decomposed  into  a  vorticity-free  field 
(Vv?  G  H-^)  plus  a  divergence-free  field  (curl  rj}  6  H).  More  generally,  if  IT  denotes  a 
projection  of  V  into  H  (recall  (2.2)),  a  projection  F  :  V  — »•  exists  such  that  F  =  /  —  IT 

(see  Temam  [45]). 

Now  suppose  that  during  a  numerical  simulation,  the  approximate  velocity  field  u^’’’ 
wanders  outside  the  space  H  of  divergence-free  fields.  We  may  then  add  to  a  “correction” 
G  so  that 

E  H  ,  or  div  u’'^  +  div  =  0 

But  since  G  there  exists  a  scalar  potential  p  such  that  =  "Vp.  Hence, 


Ap  =  —div 


(4.19) 


where  A  =  V  •  V  is  the  Laplacian.  The  scalar  p  is  termed  a  pressure  correction,  since  a  check 
with  the  projection  F  of  the  discretized  momentum  equations  into  will  show  that  'Vp  is 
related  to  VP  according  to  relations  of  the  type,  Vv?  —  AtVP  =  C>(AP),  A<  being  a  time 
step.  To  (4.19)  must  be  appended  boundary  conditions,  dpfdn  =  cr,  where  a  =  cr{h,At)  is 
selected  for  consistency  with  (2.1)  and  the  order  of  approximation  used  in  its  discretization 
in  space  and  time  (set  [58]). 

In  the  present  work,  we  shall  apply  a  pressure  correction  of  the  form  (4.19)  only  when 
the  departure  of  from  H  exceeds  a  preset  tolerance  ejoL-  Thus,  if 

|[d2y  >  ^TOL  ,  K  eTh 


we  set 

-AP  =  ArUiv  u''P 


and,  then, 


^corrected  —  +  P 


75 


5  HIGH-ORDER,  SEMI-IMPLICIT  SOLVERS 

The  regularized  Navier-Stokes  equations  can  be  written  in  the  form, 


Du 

—  =  f(«)  +  VP 


(5.1) 


where 


Du 

'm 


=  ttt  +  tt  •  V«  =  the  material  time  derivative  of  the  velocity 


F{u)  =  /  +  22/V  •  JD(u)  =  diffusion  vector 


(5.2) 


We  shall  outline  several  high-order  splitting  schemes.  For  a  survey  of  splitting  methods,  see 
Marchuk  [60]. 

5.1  High— Order  Method  of  Characteristics 


The  classical  method  of  characteristics  (though  not  actually  a  splitting  scheme)  can  be 
generalized  to  produce  a  remarkably  effective  solver.  Recall  that  the  motion  of  a  material 
fluid  particle  X  is  given  by  the  map 

X  =  x{X,t) 

where  x  is  the  spatial  position  of  X  at  time  t  and  X  is  a  smooth  invertible  function  for  every 
t  G  [0,7’].  The  characteristics  of  (5.1)  are  determined  by  the  trajectories, 


(5.3) 


and  u  can  be  written  as  a  function  of  particles  X  and  time  <  in  a  material  (Lagrangian) 
description  of  the  motion: 


u  =  u{x{X  ,t),t)  =  u{x,ty, 


Du 

Dt 


du 

dt 
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Thus,  over  a  time  interval  (nAt,(n  +  l)At), 

=  u’*(Jf)+/  [F{u{X,s))  +  VP{X,s)]ds 

Jtn 

X{X,  ^n+l)  —  x(X,<„)+  f"*' uiX,  s)ds 

Jtn 

where  u'^'^^{X)  =  u{X,tn+i)  =  ^n+i)i ^n+i)-  The  algorithm  is  thus  divided  into 

three  steps: 

1°.  Determination  of  the  Characteristics. 
set 

X"  =  X  ,  =  u{X,tn+l)  =  V°(X",0) 


As  =  {tn+i  -  tn)/N 

do  k  =  1,  ^ 

Xk  ^  X*'-^ (X*-\(A:- l)As)  As 

v’’  =  u  (X*,t„+i  —  A:As) 

m 

=  {tn+1  -  kAs) 

j=l 

end  do 


where  the  functions  tpj  =  tpj{t)  are  polynomial  functions  in  t  that  interpolate  u{x^,t)  over 
past  time  steps. 

To  further  reduce  the  error  from  the  above  procedure,  we  may  introduce  the  following 
angle  check  in  the  loop: 

Angle  Check: 


If 


<  ei  =  toll,  continue 


else  As  =  -As 

2 
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2^.  Determination 

pk  =  -  k^s) 

=  U  +  Yi^kF'‘ 

k=0 

N 

w”+'(X)  -  wqF  (u"+*(A:))  =  u  (*^,  t„)  +  Y1  ^kF’’ 

Jt=l 

3°.  Pressure-Continuity  Check 

p^+i  =  -e-'div 
0K  =  \\div 

If  0K  <S2=  <0/2,  go  to  1°, 


else 


+  p"+i 

f  Vp"‘*’^V(p  dx  =  —  f  div  c?z 


V  <p  €  M^\K) 

go  to  5“ 
end 

Here  the  functions  ip  are  the  Q’{K)  basis  functions  described  earlier;  other  notations  are  self 
explanatory. 


5.2  Multi-Step  Methods 

Multi-step  methods,  such  as  the  Adams-Bashforth  scheme,  have  been  employed  by  several 
investigators  concerned  with  spectral  approximations  of  the  Navier-Stokes  equations  (see, 
e.g.,  Canuto,  Hussani,  Quarteroni,  and  Zang  [61],  Koroczak  and  Patera  [62],  or  Ramaswamy, 
Jue,  and  Akin  [59]).  A  two-step  scheme  in  this  family  of  methods  is  indicated  as  follows: 
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Explicit  Step.  An  iV-th  order  explicit  Adams-Bashforth  step  is  embodied  in  the 
expression, 

«"+*(*)  =  «’»(*)  +  At  X]  l^kF  +  O  (At^+*) 

k=0 

where  the  are  given  in  the  following  array: 

Order  /So  0i  ^2  Pz--- 

1  1 

2  3/2  -1/2 

3  23/12  -16/12  5/12 

4  55/24  -59/24  37/24  -9/24 

etc. 

2^ .  Implicit  Correction.  The  iV-th  order  Adams-Moulton  is  expressed  by 

iV-l 

tt"+^(a;)  =  ti"  +  At  52  iSfcF  +  O  (At^+^) 

k=0 

with 

Order  ^2  ^3 

2  1/2  1/2 

3  5/12  8/12  -1/12 

4  9/24  19/24  -5/24  1/24 

The  Adams-BcLshforth/Adams-Moulton  schemes  can  be  used  in  a  predictor-corrector  algo¬ 
rithm  with 

iv-i 

u*  =  tt"  -h  At  52  (predictor) 

k=0 

_ 

=  ti"  -|-  At^oF(u*)  -f  At  52  (corrector) 

k=l 

5.3  Runge-Kutta  Schemes 

The  cljissical  explicit  Runge-Kutta  schemes  provide  another  approach  toward  the  develop¬ 
ment  of  high-order  flow  solvers.  In  general,  these  advance  the  solution  in  time  according 
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=  u"'  +  At^2 


£  =  1 


ki  =  F  +  At  ^2  ^ijkjj 

with  1  <  m  <  i  —  1  for  explicit  schemes,  m  =  i  for  semi-implicit  schemes,  and  i  +  l  <  m  <  N 
for  fully  implicit  schemes.  See  Butcher  [63]  for  full  details  and  tables  of  the  coefficients  Ujj. 

We  compare  results  obtained  with  some  of  these  schemes  in  Section  7. 


6  ERROR  ESTIMATION  AND  ADAPTIVITY 

The  /ip-adaptive  strategy  must  involve  a  running  audit  of  the  error  in  the  approximate 
velocity  and  pressure  fields.  Existing  mathematical  theories  for  a  posteriori  error  estimation 
primarily  apply  to  simple  elliptic,  scalar  boundary-value  problems  and  few  results  exist  for 
error  estimation  of  solutions  to  the  Navier-Stokes  equations.  A  study  of  error  estimation 
techniques  for  /ip-finite  element  approximations  of  elliptic  problems  is  given  in  [42];  some 
extensions  of  these  ideas  to  time-dependent  problems  are  discussed  in  [63]  and  error  residual 
methods  are  adapted  to  calculations  of  compressible  flow  problems  in  [64]. 

In  this  investigation,  we  extend  the  error  residual  methods  of  Ainsworth  and  Oden  [65] 
to  incompressible  Navier-Stokes  equations  by  exploiting  the  particular  structure  of  the  ap¬ 
proximations  inherent  in  the  splitting  algorithms  described  in  the  preceding  section. 

6.1  A  Posteriori  Error  Estimators 


The  (n  -f  l)st-step  of  a  semi-explicit  method  of  the  type  discussed  in  Section  5  requires  the 
solution  of  a  linear  elliptic  system  of  the  form, 

A<e-iV(V-u"+^) 

(6-1) 


=  Atr^^+u 


"+1  -L  «»»+7 


with,  for  example, 


=  u”  (x^) 

for  the  method  of  characteristics.  Setting 

=  u  ;  AD{u)  =  At  (2i/D{u)  +  e“'l  V  •  u) 
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the  problem  reduces  to  the  symmetric,  linear,  elliptic  boundary-value  problem, 


—V  •  AD{u)  -1-  u 

=  f 

in 

n 

u 

=  0 

on 

Ti 

n  •  Au 

=  9 

on 

r2 

(6.2) 


Here  A  is  the  fourth  order  tensor,  Ajkt  =  At{2i/SnSjk  +  2vSjk^jt  +  £  The  corre¬ 

sponding  bilinear  and  linear  forms  are 


a :  y  X  y 


,  L:  V^R 


2(tt,  r)  =  /  (^D{v)^ AD{u)  +  u  •  dx  =  ^  akc{u,v) 

Ken 

(6.3) 

L{v)  =  f  f  -v  dx-\-  f  g  V  ds  =  ^  Lk{v) 

JU  JTi 

(6.4) 

Here  d^i-  ,  •)  and  L}({-)  are  the  bilinear  and  linear  forms  defined  on  a 
(recall  (2.6),  (2.7)).  We  define  the  local  and  global  energy  norms  by 

partition  Th  of  0 

11«IIe  =  \/5(w,w)  ,  ||M||£;.if  =  yJdK{u,u) 

(6.5) 

Thus 

We  wish  to  estimate  the  approximation  error  of  e(®,<)  =  u{x,t)  —  u^^(x,t)  a.t  t  = 
(n  -{-  1)A<  to  within  term  of  in  the  energy  norms  (6.5). 

We  next  introduce  the  following  notations  and  conventions  (Cf  [65]); 

_  N  _ 

•  The  partition  of  fl  has  N  =  N{h)  elements  K  =  CIk,  H  =  UfA'  G  7/i}  =  (J 

K=l 

(thus  “K”  is  used  to  indicate  both  a  subdomain  and  an  index  of  the  subdomain) 
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•  Fa'i,  =  0  <  K,L  <  N  are  sets  consisting  of  a  finite  number  p{I\,L)  of 

segments  such  that 


^KL 


p{K,L) 

u 


M=l 


Sa'L  =  isolated  (nodal)  points  (Sa'l  P)  “  ®) 


•  Tok  =  dnKf]dn ,  5^;,  =  U 

L,M 


•  E  =  E  {Th)  =boundary  segments  of  the  partition 

=  U  r^L 

K.LsO 

K>L 

l<M<p(K,L) 


•  El  =  segments  lying  interior  to  Q 


•  <^KL  =  —ctlk 


r  +1  if  K>L 
{  -1  if  K  <L 


N 


u 

K,L=1 

K>L 


pA/ 

^  KL 


•  n{s)  =  <TKLnK{s)  =  (TLKniis)  ,  S  e 

thus  n  is  a  unit  normal  pointing  outward  from  the  subdomain  with  largest  index 


•  ^KL  •  ^KL  i  <^Kl(^)  +  —  I5 

0<K,L<N,  1<M<  p(K,L) 


•  Jumps  and  convex  combinations: 


[r  •  n|pM  =  ctklI'^k  ■ 

i\  L 

(*’  •  '^)air«^  =  0‘KLl'f’K  ■  riK  -  QlkJ^l  •  riL 

(')Vii  —  trace  of  v  ^  onto  dn^,  [r  •  n|  is  the  jump  in  v  •  n  across  neighboring 

subdomains,  while  {v  •  n)a  is  a  linear  weighted  average  of  v  •  n  between  neighboring 
subdomains) . 
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The  element  residual  functions  are  then  defined  by, 

ri<{x)  =  fix)  +[V-AD  (u'‘p)  -  u'‘p)  (x)  ,  x  e  K 

g  -Tifc  •  AD  s  e  r^fon^2  ,  1  <  <  /j(A',0) 

-aKL  [[n  .  AD  (u'*p)]]  ,  s  e  ,  I  <  K,  L  <  N 

1  <  M  <  piK,L) 

RKds 

Ok  =  1 

[  0  if  otherwise 

Each  finite  element  is  now  endowed  with  a  complimentary  energy  given  by  the  functional 
GkCo-)  =  tr^A~^<T  dx  -  ^  ^  I V  •  <r  +  tk\^  dx  (6.7) 

Theorem  6.1  (Cf  [60]).  Let  e  denote  the  approximation  error 

e  —  u  —  (6-S) 

where  u  is  the  exact  solution  to  (6.1)  and  is  its  hp-finite  element  approximation.  Then 

Ml  <  -2  E  dKior)  (6.9) 

KO-h 

for  all  (T  G  VE/c,  where 

=  {<7i,(x)G(L2(/0)"'|V-<re(A2(/v))", 

cr.j  =  <Tj,  ,  1  <  i,j  <  3  ,  (TijUKi  =  Rhiis),  (6.10) 

on  5/v  \ri} 


f  meas{K)  ^  rKdx  +  ^ 


I  -ddKPidn  =  0 


83 


Defining, 


Vi'ifr)  =  4(<r)  + 

e\'{(r)  =  f  A~^cr  dx 

J  hC 

=  /  \V  ■<T-^rK\^dx 

**  K  j 


(6.11) 


we  see 


that 


Mil  <  Y1  ^Ki^)  n  (6-12) 

K  eTh  K  eT), 

Thus,  the  functionals  r]f^(<r)  provide  a  bound  on  the  global  error  for  any  “pseudo  stress”  <r. 
It  remains  only  to  construct  an  efficient  way  to  calculate  (or  approximate)  a  suitable  local 
field  er. 

One  technique  that  leads  to  asymptotically  exact  estimates  is  as  follows; 

Let  U  €  be  such  that 


/  D{x)'^ AD{U)dx  =  f  {rK  +  6K)-xdx 
*'  /C  J  K 

+  /  RK-Xds+  [  g-xds 

JdK\ri  JdKnr2 

yx^V^'P+i{K) 

It  can  be  shown  that  there  exist  a  unique  U  satisfying  this  equation.  Set 


(6.13) 


=  AD{U) 


(6.14) 


Then  it  can  be  shown  that  A^-(<r*)  <  A^-(<t),  V  <t  €  Wk,  and  that  is,  therefore,,  an 

acceptable  error  indicator.  We  choose 


^‘dx  +  meas{K)6j^ 


(6.15) 
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6.2  Adaptive  Strategy 

The  adaptive  algorithm  employed  is  briefly  outlined  is  as  follows: 


1.  Begin  with  an  initial  mesh,  generally  obtained  after  solving  the  Stokes  problem  with 
one  or  two  /i-refinements  and  uniform  p,  or  with  p  higher  in  regions  where  effects  of 
singularities  are  expected  to  be  significant. 

2.  Solve  the  problem  on  the  initial  mesh. 

3.  Compute  due  to  the  increase  of  every  unconstrained  degree  of  freedom  of  element 
K. 

4.  Choose  the  change  in  or  that  produced  the  maximum  change  in  error 

At]^  =  max 

and  redistribute  hx  and  accordingly. 

5.  Go  to  item  2. 

Details  of  this  procedure  are  given  in  [43]. 


7  NUMERICAL  EXPERIMENTS 

In  this  section,  we  present  selected  results  of  several  numerical  experiments  designed  to 
test  the  methods  and  the  general  approach  described  earlier.  These  primarily  concern  stan¬ 
dard  model  problems,  such  as  the  backstep  channel  problem  and  the  driven  cavity  problem. 

7.1  The  Backstep  Channel  Problem  in  Two  Dimensions 

A  standard  benchmark  problem  is  the  so-called  backstep  channel  problem  depicted  in  Fig. 
5.  The  inflow  velocity  is  parabolic, 


and  no-flow  conditions  are  imposed  at  all  other  boundaries  except  the  outflow  boundary 
xo  =  &  -f  c,  where  we  set  (r{u{t),  P{t))  •  n  =  0.  The  recirculation  lengths  shown 

depend  upon  the  Reynolds  number,  and  data  on  the  variation  of  Li  and  I2  with  Re  has 
been  obtained  by  several  investigators  using  other  schemes  (see,  in  particular,  Ghia.  et  al. 
[64]). 
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Figure  5:  Geometry  and  notation  for  the  two-dimensional  backstop  channel  problem. 
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The  lirst  results  cited  here  focus  on  our  comparisons  of  the  schemes  described  in  Section 
5.  No  attempt  was  made  to  optimize  or  vectorize  the  coding  in  these  calculations,  and  a 
GMRES  scheme  was  employed  to  solve  the  diffusion  steps  in  all  cases.  The  uniform  mesh 
of  IQ^  elements  shown  in  Fig.  6,  corresponding  to  around  1200  degrees  of  freedom,  was 
used. 

Computed  velocity  profiles  over  a  portion  of  the  mesh  are  shown  in  Figs.  7  and  8.  Figure 
7  contains  results  obtained  using  a  fourth-order  Runge-Kutta  (RK4)  scheme.  There  we 
see  computed  velocity  fields  at  t  =  2.00  secs,  obtained  using  a  time  step  At  =  0.1  sec. 
Figure  8  contains  similar  results  for  the  same  time  step  At  but  computed  using  a  first-order 
MOC  (method  of  characteristics)  technique.  The  computed  velocity  fields  are  essentially  the 
same,  but  were  obtained  by  the  MOC  required  only  40  percent  of  the  CPU  time  of  the  RK4 
technique. 

Similar  results  for  a.  Re  =  600  are  shown  in  Figs.  9-11.  In  Fig.  9,  we  see  computed 
velocities  delivered  by  the  RK4  scheme  at  t  =  1.00  secs,  corresponding  to  At  =  0.1  sec.  A 
high-order  (fourth-order)  MOC  scheme  was  also  used  to  solve  this  problem  and  results  are 
shown  in  Fig.  10.  This  calculation  required  approximately  75  percent  of  the  CPU  time  used 
by  the  RK4  scheme.  A  fourth-order  Adams-Bashforth  scheme  was  also  used  to  solve  this 
example,  but  it  failed  for  At  =  0.1  sec.  and  was  stable  for  At  <  0.05  sec.  Results  obtained 
for  this  scheme  at  t  =  0.50  sec.  are  shown  in  Fig.  11. 

These  calculations  were  repeated  for  Re  =  100  and  1000,  and  a  comparison  of  results  is 
given  in  Table  1  below.  Variances  in  time  steps  shown  in  the  table  are  due  to  the  conditional 
stability  of  -the  explicit  steps  in  the  RK4  scheme  and  the  fourth  order  Adams- Bashorth 
scheme,  which  proved  to  be  unstable  at  At  =  0.05  sec.  for  Re  =  1000.  The  method  of 
characteristics  appears  to  be  unconditionally  stable  and  to  be  superior  to  the  other  high- 
order  schemes  for  these  calculations.  The  relative  efficiency  of  these  schemes  could,  of  course, 
change  with  appropriate  optimization  of  the  coding.  All  calculations  were  run  on  a  single 
processor  of  the  Alliant  FX8  computer  and  are  normalized  with  respect  to  the  MOC  for 
Re  =  100.  Results  obtained  using  multiple  processors  are  underway  and  will  be  reported  in 
future  work. 


Table  1:  Comparison  of  High— Order  Schemes 
Q^/Q^  Elements  1200  degrees  of  freedom 


Re 

300 

600 

At 

CPU 

At 

CPU 

RK4 

0.1 

2.43 

0.1 

2.44 

Adams/B4 

0.1 

3.02 

0.05 

— 

MOC4 

0.1 

1.00 

0.1 

1.04 
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Figure  6:  A  uniform  mesh  of  Q^IQ^  elements  with  1202  degrees  of  freedom;  a/h 
cja  =  1/24,  Re  —  300  and  Re  =  600. 
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Figure  7;  Velocity  profiles:  fourth-order  Runge-Kutta  scheme,  elements, 

t  =  2.0  sec,  Re  =  300. 
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Figure  10:  Computed  velocity  fields  at  t  =  1.00  sec,  fourth-order  MOC,  At  =  0.1  sec, 
Re  =  600. 
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Figure  11;  Computed  velocity  field  at  f  =  0.50  sec;  conditionally  stable,  fourth-order 
Adams- Bashforth  scheme  with  At  =  0.05  sec,  Re  =  600. 
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Q'^IQ^  Elements  1756  degrees  of  freedom 


Re 

100 

1000 

mmm 

At 

CPU 

At 

CPU 

RK4 

0.05 

6.60 

0.01 

6.15 

Adams/B4 

0.05 

5.66 

0.01 

mi 

MOC4 

0.1 

1.00 

0.1 

2.57 

Table  2  contains  computed  values  of  the  reattachment  length  Li  and  L2  of  Fig.  5  com¬ 
puted  using  a  fixed  grid  of  Q^fQ^  elements  for  various  Reynolds  numbers  for  a  =  0.51485 
compared  with  those  obtained  by  Ghia,  et  al.  [64]  using  a  12,870  degrees  of  freedom  finite 
difference  calculation.  These  results  correspond  to  solutions  at  t  =  100.0  sec.  and  were 
obtained  using  only  1394  degrees  of  freedom. 

Table  2:  Reattachment  Length  Comparisons 


Re 

Li 

Li 

Present 

Ref.  [64] 

Present 

Ref.  [64] 

150 

3.23 

3.38 

— 

— 

300 

4.96 

4.95 

4.05 

3.75 

507 

6.08 

6.25 

4.70 

4.75 

600 

6.39 

6.50 

4.97 

5.00 

7.2  A  Backstep  Channel  in  Three  Dimensions 

Some  preliminary  results  on  a  three-dimensional  backstep  channel  calculation  are  displayed 
in  Fig.  12  which  contains  ui-velocity  profiles  for  the  Stokes  problem,  the  solution  of  which 
is  used  as  the  initial  conditions  for  integrating  the  Navier-Stokes  equations.  Here  a  uniform 
mesh  of  Q^IQ^  elements  is  used  (it  being  understood  that  Q''  corresponds  to  tensor  products 
of  polynomials  of  degree  r  on  a  hexahedral  element).  Solutions  for  Re  =  300,  40  elements, 
At  =  0.1  sec.  at  t  =  5.0  sec.  are  shown.  A  more  detailed  analysis  of  this  problem  is  underway 
and  is  to  be  discussed  in  a  future  report. 


7.3  Flow  Around  a  Square  Body,  Re  =  250 

Figure  13  shows  an  initial  hp  mesh  of  cubic  and  quintic  elements  surround  a  square  obstacle 
in  a  rectangular  flow  channel.  A  parabolic  inflow  velocity  is  prescribed,  the  intensity  of  which 
corresponds  to  a  Reynolds  number  of  250  relative  to  the  side  length  of  the  square.  The  finite 
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Figure  12:  Computed  velocity  field  in  a  three-dimensional  backstop  channel:  Re  —  300, 
t  =  5.0  sec.,  MOCl. 
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element  model  contains  only  1656  degrees  of  freedom.  A  second-order  semi-implicit  Runge- 
Kutta  scheme  was  used  to  integrate  the  equations  of  motion,  with  penalty  parameter  £  set 
uniformly  to  lO-"*. 

Computed  velocity  vectors  and  contours  of  the  longitudinal  coinponent  of  velocity  are 
shown  in  Fig.  14  at  <  =  230  sec.  Contours  of  the  U2-velocity  at  this  time  are  shown  in  Fig. 
15. 


7.4  Two-Dimensional  Driven  Cavity,  Re  =  5000 

The  classical  driven  cavity  problem  is  next  considered  as  a  standard  benchmark  test.  Two 
cases  are  presented.  In  the  first  case,  only  ^-refinement  is  used  with  p  =  4  uniformly  over 
the  square  grid.  The  mesh  is  shown  in  Fig.  16,  has  2393  degrees  of  freedom.  The  Ui-velocity 
of  the  “lid”  is  prescribed  to  be  uniform  and  is  set  to  correspond  to  a  Reynolds  number  of 
5000.  The  adaptive  algorithm  automatically  refines  the  mesh  near  the  corners  to  resolve  the 
corner  singularities. 

In  this  test,  the  method  of  characteristics  is  used  to  integrate  the  solution  to  an  apparent 
steady-state  solution.  Computed  ui  =  u  and  vi  =  v  velocity  contours  are  shown  in  Fig. 
17.  Computed  pressure  contours  are  shown  in  Fig.  18  compared  with  those  of  Gresho  and 
Chan  [58],  which  were  obtained  using  a  fine  graded  mesh  of  bilinear  elements.  It  is  observed 
that  excellent  agreement  is  achieved.  Notice  that  the  plotting  routine  has  detected  the  slight 
discontinuities  in  the  pressure  field,  indicated  by  small  but  stable  and  non-oscillatory  jumps 
at  boundaries  of  the  larger  elements. 

A  second  case  was  run  using  the  nonuniform  hp  mesh  shown  in  Fig.  19  which  involved 
only  1442  degrees  of  freedom.  In  this  case,  the  adaptive  procedure  produced  the  mesh  shown 
with  Q^IQ^,  and  Q^IQ^  elements.  Again,  the  MOCl  scheme  was  used  for  Re  =  5000. 

Results  essentially  coincided  with  those  in  Figs.  17  and  18.  To  exaggerated  discontinuities 
in  the  pressure,  a  coarser  contour  gridding  was  used  and  the  resulting  pressure  contours  are 
shown  in  Fig.  20. 

7.5  Three-Dimensional  Driven  Cavity 

A  three-dimensional  version  of  the  driven  cavity  problem  was  also  tested  and  results  are 
shown  in  Fig.  21. 

A  relatively  crude  mesh  of  22  cubic  elements  (Q^/Q^)  w'th  one  level  of  /i-refinement  is 
shown  in  Fig.  21.  This  corresponded  to  a  model  with  only  2244  degrees  of  freedom.  A 
uniform  velocity  U2  =  1.0  is  presented  over  the  lid  face  (a:  =  0)  and  the  RK4  scheme  is  used 
to  integrate  the  Navier-Stokes  equations.  In  this  case.  Re  =  100  and  =  0.1  sec.  The 


96 


Figure  13:  hp  mesh  for  simulation  of  flow  past  a  square  obstacle;  Re  —  250. 
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Figure  16;  An  hp  mesh  for  the  classical  driven  cavity  problem,  with  adapted  /i-refinements 
and  p  =  4. 
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computed  profiles  at  <  =  10.0  At  are  shown  superimposed  on  the  mesh.  A  detailed  discussion 
of  results  obtained  on  this  and  other  three-dimensional  benchmark  problems  is  forthcoming. 
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